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FOREWORD 


This report is a revision of the document “MASTERFIT” — 1986, dated August 1, 1986, which it 
supersedes. The model changes during 1986 and 1987 were rather minor. They included small correc- 
tions for feed rotation and refraction in modeling antenna axis offsets (Section 2.11, both previously 
neglected). Section 2.7 now includes the contributions from the dependence on the UT 1 transforma- 
tion to the nutation partials. A temporary optional model to correct semiannual and annual nutation 
terms has been installed. Finally, it is now possible to estimate the surface temperature when using 
the CfA mapping function (Section 4.3). Figures 1 and 2 have been redrawn to be more realistic 
representations of the geometric delay. A number of misprints have also been corrected. The present 
document corresponds to MASTERFIT version 358, which has been in use since September, 1987. 
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ABSTRACT 


This report is a revision of the document of the same title (1986), dated August 1, 1986, which it 
supersedes. Model changes during 1986 and 1987 included corrections for antenna feed rotation, re- 
fraction in modeling antenna axis offsets, and an option to employ improved values of the semiannual 
and annual nutation amplitudes. Partial derivatives of the observables with respect to an additional 
parameter (surface temperature) are now available. New versions of two figures representing the geo- 
metric delay are incorporated. The expressions for the partial derivatives with respect to the nutation 
parameters have been corrected to include contributions from the dependence of UT1 on nutation. 
The authors hope to publish revisions of this document in the future, as modeling improvements 
warrant. 
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SECTION 1 


INTRODUCTION 

In applications of radio interferometry to geodynamics and astrometry, observed values of delay 
and delay rate obtained from observations of many different radio sources must be passed simultane- 
ously through a multiparameter estimation routine to extract the significant model parameters. As 
the accuracy of radio interferometry has improved, increasingly complete models for the delay and 
delay rate observables have been developed. This report describes the current status of the delay 
model used in the Jet Propulsion Laboratory multiparameter estimation program “M ASTERFIT* . It 
is assumed that the reader has at least a cursory knowledge of the principles of VLBI. 

The delay model is the sum of four major model components: geometry, clock, troposphere, 
and ionosphere. Sections 2 through 5 present our current models for these components, as well as 
their partial derivatives with respect to parameters that are to be adjusted by multiparameter fits 
to the data. The longest section (2) deals with the geometric delay and covers the topics of time 
definitions, tidal effects, coordinate frames, Earth orientation (universal time and polar motion), 
nutation, precession, Earth orbital motion, wave front curvature, gravitational bending, and antenna 
offsets. Section 6 describes the technique used to obtain the delay rate model from the delay model. 
Section 7 gives the values of physical constants used in MASTERFIT, while section 8 outlines model 
improvements that may be required by more accurate data in the future. 
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SECTION 2 


GEOMETRIC DELAY 


The geometric delay is that interferometer delay which would be measured by perfect instrumen- 
tation, perfectly synchronised, if there were a perfect vacuum between the sources observed and the 
instrumentation. Typically this delay can be as great as 0.02 second for Earth-fixed baselines, and 
changes rapidly (up to 1.5 /xsec per second) as the Earth rotates. In general this component of the 
delay is by far the largest component of the observed delay. The main complexity of this portion 
of the model arises from the numerous coordinate transformations necessary to relate the reference 
frame used for locating the radio sources to the Earth-fixed reference frame in which station locations 
are represented. 

In the following we will assume, unless otherwise stated, that “celestial reference frame” means 
a reference frame in which there is no net proper motion of the extragalactic radio objects which 
are observed by the interferometer. This is only an approximation to some truly “inertial* frame. 
Currently, this celestial frame implies a geocentric, equatorial frame with the equator and equinox of 
J2000 as defined by the 1976 IAU conventions with the 1980 nutation series (Seidelmann, 1982, and 
Kaplan, 1981). 

In this equatorial frame, some definition of the origin of right ascension must be made. We will 
not discuss that in this paper, since one definition is at most a rotation from some other definition, 
and can be applied at any time. The important point is that consistent definitions must be used 
throughout the model development. The need for this consistency will, in all probability, eventually 
lead to our defining the origin of right ascension by means of the JPL planetary ephemerides, followed 
by our using interferometric observations of both natural radio sources and spacecraft at planetary 
encounters as a means of connecting the planetary and the radio reference frames. 

Also, unless otherwise stated, we will mean by “terrestrial reference frame” some reference frame 
tied to the mean surface features of the Earth. Currently, we are using a right-handed version of 
the CIO reference system with the pole defined by the 1903.0 pole. In practice, this is accomplished 
by defining the position of one of the interferometric observing stations (generally DSS 14 at the 
Goldstone Deep Space tracking complex), and then by measuring the positions of the other stations 
under a constraint. This constraint is that the determinations of Earth orientation agree on the average 
with the Bureau International de l’Heure (BIH) (1982) measurements of the Earth’s orientation over 
some substantial time interval (« years). This procedure, or its functional equivalent, is necessary 
since the interferometer is sensitive only to the baseline vector as measured in the celestial frame. The 
VLBI technique is purely geometric and does not have any preferred origin point directly based on 
the structure of the Earth. The rotation of the Earth does, however, provide a preferred direction in 
space which can be associated indirectly with the surface features of the Earth. 

In contrast, geodetic techniques which involve the use of satellites, or the Moon, are sensitive 
to the center of mass of the Earth as well as the spin axis. Thus, those techniques require only a 
definition of the origin of longitude. We anticipate that laser ranging to the retroreflectors on the 
Moon (LLR) will allow a realizable practical definition of a terrestrial frame, accurately positioned 
relative to a celestial frame tied to the planetary ephemerides. Co-location of the laser stations with 
the interferometer elements will be required to achieve this, however. 

Careful definitions and experiments of this sort will be required to realize a coordinate system of 
centimeter accuracy. In the meantime, we must establish interim coordinate systems carefully enough 
so that we do not degrade the intrinsic accuracy of the interferometer data by introducing “model 
noise” . 

Several equivalent formulations are possible for the calculation of geometric delay. Our interest 
in also using this technique for radio sources within the Solar System has led to the formulation 
presented here. The problem of calculating the time delay for these “local” sources is simplified if the 
calculations are carried out in a frame at rest relative to the center of mass of the Solar System. For 
our present data, Lorentz transformations adequately account for the relativistic effects which arise 
from the relative motion of the objects involved. General relativity is required only for those effects 
which are associated with the bending of ray paths by the more massive objects in the Solar System. 
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For the above reasons we have performed the calculation as outlined below: 

1. Specify the coordinates of the two stations as measured in an Earth-fixed frame at the time that 
the wave front intersects station # 1. Let this time be t\ as measured by a clock in the Earth-fixed 
frame. 

2. Modify the station locations for Earth-fixed effects such as solid Earth tides, tectonic motion, 
and other local station motion. 

3. Transform these station locations to a celestial coordinate system with the origin at the center of 
the Earth, but moving with the Earth. This is a composite of 10 separate rotations, represented 
by a rotation matrix Q(t). 

4. Perform a Lorentz transformation of these station locations from the Earth-fixed frame to a 
frame at rest relative to the center of mass of the Solar System, and rotationally aligned with the 
celestial geocentric frame. 

5. In this Solar System barycentric frame, compute the time interval for the passage of the specified 
wave front from station #1 to station #2. Include the effects of the curvature of the wave front 
due to the finite distance of the source if relevant. Also, add in the effects of the differential 
gravitational retardation of the signal received by way of the two separate interferometer ray 
paths. 

6. Perform a Lorentz transformation of this calculated geometric delay back to the celestial geocen- 
tric frame moving with the Earth. Except for negligible effects discussed later, this produces the 
model for the geometric portion of the observed delay. 

7. To this geometric delay, add the contributions due to clock offsets, to tropospheric delays, and 
to the effects of the ionosphere on the signal (see sections 3 through 5 below). 

The remainder of this section provides the details for the first six steps of the general outline 
above. Note that we have invoked only special relativistic coordinate transformations. Gravitational 
effects are incorporated as correction terms added on to the geometric delay once it is calculated in 
the Solar System barycenter. 


2.1 TIME INTERVAL FOR THE PASSAGE OF A WAVE FRONT 
BETWEEN TWO STATIONS 


The fundamental part of the geometric model is the calculation (step #5 above) of the time 
interval for the passage of a wave front from station #1 to station #2. We actually do that calculation 
in a coordinate frame at rest relative to the center of mass of the Solar System. However, the basic 
formulation is independent of the frame used. This part of the model is presented first to provide a 
context for the subsequent sections, all of which are heavily involved with the details of time definitions 
and coordinate transformations. We will use the same subscript and superscript notation which is 
used in section 2.10 below to refer to the station locations as seen by an observer at rest relative to 
the center of mass of the Solar System. 

First, we calculate the interval that would be observed if the wave front were planar. Next, we 
generalise this calculation to a curved wave front, and finally, we take into account the incremental 
effect which results from the fact that we must consider wave fronts that propagate through the various 
gravitational potential wells in the Solar System. To start, consider the case of a plane wave moving 
in the direction, k, with station 2 having a mean velocity, fi 2 , as shown in figure 1. The time delay 
is the time it takes the wave front to move the distance l. Measured in units of time, this distance is 
the sum of the two solid lines perpendicular to the wave front in figure 1: 

t; - h = k • [r 2 (t!) - ti (ti)] + k • fi 2 [t* - h] (2.1) 


This leads to the following expression for the geometric delay: 


= £-M«i)-*i(*i)1 

1-k fi 2 


( 2 . 2 ) 
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In the case of a signal generated by a radio source within the Solar System it is necessary to 
include the effect of the curvature of the wave front. As depicted in figure 2, let a source irradiate two 
Earth-fixed stations whose positions are given by r<(t) relative to the Earth’s center. The position of 
the Earth’s center, R c (ti), as a function of signal reception time, 1 1 , at station #1 is measured relative 
to the position of the emitter at the time, t e , of emission of the signal received at time tj.. While this 
calculation is actually done in the Solar System barycentric coordinate system, the development that 
follows is by no means restricted in applicability to that frame. 

Suppose that a wave front which leaves the source at time t e reaches station #1 at time fi and 
arrives at station #2 at time t\. The geometric delay in this frame will be given by: 

r = t5-*i = |Ba(*5)|-|lli(ti)| (2.3) 



IT AT TIME t, 


Figure 1. Geometry for calculating the transit time of a plane wave front 


where all distances are 

again measured in units of light travel time. If we approximate the velocity of 

station #2 by 

a ^-2 (* 2 ) — H- 2 (*l) 

h ~ *5 — 

(2.4) 

and use the relation 

Rn(ti) = R c (ti) + r»(ti) 

(2.5) 

we obtain: 

r = |R c (ti) + ra(ti) + fi 2 r \ ~ PM*i) + r i(*i)| 



= Rc{t 1 ) [ |Rc + «a( — |Rc + *i| ] 

(2.6) 


4 



where 


(2.7) 


and 


^2(^1) + j 9 2 T 
R c (t 1 ) 


«i 


£lM 

Jle(tl) 


( 2 . 8 ) 


For ei and c 2 < 10 -4 , we need to keep only terms of order e 3 in a sixteen-place machine in order to 
expand the expression for r in equation (2.6). This gives us: 


R c fofri) ~ri(<i)l + Rc^c(r) 
(l-Rc-^a) 2 [1 — R c ^ 2 ] 


(2.9) 



Figure 2. Geometry for calculating the transit time of a curved wave front 


where to order e 3 

A e (r) = [e% - ej] - [(R c e 2 ) 2 + (Re *i) 2 + (R c « 2 ) 3 ~ (R c » 2 )e 2 ~ (R c «i) 3 + (R c «i)ef] (2.10) 

The first term in (2.9) is just the plane wave approximation, t.e., as R c — ► 00 , R c — ♦ k, with the 
second term in brackets in (2.10) approaching zero as r 2 /R c . Given that the ratio of the first term 
to the second term is w r/i? c , wave front curvature is not calculable in a sixteen-place machine for 
R > 10 16 X r. For Earth-fixed baselines, the requirement that the effects of curvature be less than 0.01 
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cm implies that the above formulation (2.10) must be used for R < 3.6 x 10 14 km, or approximately 
2 X 10 6 astronomical units. 

The procedure for the solution of (2.9) is iterative for e < 10 -4 , using the following: 


r n 


R c A e (r n - X ) 

T ° 2[1> — R/9 2 ] 


where 


TO — Tpi ane mave 


For e > 10 ~ 4 , directly iterate on the equation (2.6) itself, using the procedure: 


( 2 . 11 ) 


( 2 . 12 ) 


Tn = R«|H.c + *2(fn-l)| ~ Rc|Re + *1 1 


(2.13) 


where again to is the plane wave approximation. 

Because of the retardation of a light signal propagating in a gravitational potential well, the 
computed value for the differential time of arrival of the signals at ri(ti) and r 2 (t 2 ) must be incre- 
mented from the value calculated above. For the geometry illustrated in figure 3, the required delay 
increment, A < 3 P , is given by Moyer (1971): 


A Gp = 


(1 + 1ppn)p p 


t, + r 2 (t 2 ) + r,a 
r, + r 2 (t$) - r, 2 


— In 


r» + r i(tj) + r,i 
r. + ri(ti) - r,i 



(2.14) 


where r,i is defined as: 


r »i — |l"»(t«) l"«(^e)| 


(2.15) 


As will be shown below, some care must be taken in defining the positions given by r,, r 2 (t5), and 
n(ti). We have chosen as the origin the position of the gravitational mass at the time of closest 
approach of the received signal to that object. The position, r,, of the source relative to this origin is 
the position of that source at the time, t e , of the emission of the received signal. Likewise, the position, 
r«(t»), of the i th receiver is its position in this coordinate system at the time of reception of the signal. 
Even with this care in the definition of the relative positions, we are making an approximation, and 
implicitly assuming that such an approximation is no worse than the approximations used by Moyer 
(1971) to obtain (2.14). Here ippn is the 7 factor in the parameterized post-Newtonian gravitational 
theory: 

lPPN = (2 - 16) 

where w is the coupling constant of the scalar field. For general relativity, 'Ippn = 1, t.e., w — ► 00 . 
However, we allow 7 ppff to be a “solve-for” constant so that by setting 'Ippn — —1, also have 
the option of “turning off” the effects of general relativity on the estimate of the delay. This proves 
useful for software development. The gravitational constant, /i p , is 


H P = Grrip 


(2.17) 


where G is the universal gravitational constant, and m p is the mass of the body. Dropping the time 
arguments, we have: 


A Gp = 


(1 + 1ppn)^p 


■ In 




r, + r 2 + r, 2 
r, + ri + r 4l 


r, + ri - r tl 
r, + r 2 - r .2 


(2.18) 
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This formulation is fine for r, » u » r,i, but can be put in a computationally better form for 
the case of distant sources with closely spaced VLBI receivers, ».e., |r 2 — ri |/ri — ► 0 , r</r, —* 0. For 
these sources, expand A g p in terms of n/r,, r,i/r,, and make use of the relationship 

r.i = [ r ] - 2r, ■ n + r?] 1/2 w r 4 [l - r< • ?,] 

This leads to 


Acp = 


(1 + 7ppn)Pp 


In 


ri + ri • r, 

f2 + r 2 ■ ? , 


for rijr, -* 0. 


(2.19) 

( 2 . 20 ) 



Figure 3. A schematic representation of the geodesic connecting two points in the 
presence of a gravitational mass 


If we further require that jr 2 — ril/jr^ — » 0, and make use of 

r 2 = ri + Ar 

then: 

r 2 +r 2 •?, = ri [l + 2 • Ar/r! + (Ar/ri) 2 j ^ + ri-r, +Ar-r, 

» (l + ?i ■ Ar/rij + r x • r, + Ar • r. 

In the limit of Ar/ ri — + 0: 

r 2 [l + r 2 • r,j -»ri[l + ri r,] + Ar ■ [?i + r.] 
Substituting into (2.20) and expanding the logarithm, we obtain: 

A (1 + '1ppn)Pp [^2 - rj] ■ [fi +r,] 

° p_ ' ri[l + fi • r,] 


( 2 . 21 ) 

( 2 . 22 ) 

(2.23) 

(2.24) 
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Using whichever of these three formulations (2.18, 2.20 or 2.24) is computationally appropriate, 
the model calculates Ac P for each of the major bodies in the Solar System (Sun, planets, Earth, and 
Moon). The total gravitational correction used is: 


N 

Ac = A Gp 

p= 1 


(2.25) 


VELOCITY 


EMITTER 



Figure 4 . A schematic representation of the motion of a gravitating object during the 
transit time of a signal from the point of closest approach to reception by 
an antenna 


where the summation to N is over the major bodies in the Solar System. As mentioned above, 
care must be taken to use the positions of the emitter, the gravitational object, and the receivers at 
the appropriate times. For a grazing ray emitted by a source at infinity, using the position of the 
gravitating body G at the time of reception of the signal at station #1 rather than at the time of closest 
approach of the signal to G can cause a 15-cm error on baselines with a length of one Earth radius as 
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shown by the following calculation. Prom figure 4, the calculated distance of closest approach, R, 
changes during the light transit time, tught transit, of a signal from a gravitational object at a distance 
Reg by: 


AJR W Reg © * bight transit — © * R%g/ ^ 


(2.26) 

Since the deflection (see Appendix A) is: 





(2.27) 

MA6) = -A6[ A /] = Ae[®*- ] - Ae[*“] 

ISO 
1 dt 

(2.28) 


We consider the two bodies of largest mass in the Solar System: the Sun and Jupiter. Their respective 
deflections A0 are 8480 and 73 nanoradians. The angular velocities are estimated to be 6 x 10 -11 
and 1.7 X 10 -8 rad/sec for the Sun and Jupiter. Note that Eq. (2.26) does not apply to the Sun. 
Its motion in the barycentric frame has a period of 11 years with a radius of the order of the Sun’s 
radius. Using approximate radii and distances from Earth to estimate Reg and 6, Eq. (2.28) gives 
2.5 X 10 -8 rad for Jupiter; the corresponding value for the Sun is 7 x 10 -11 . For a baseline whose 
length equals the radius of the Earth, 8(A0)Re is thus approximately 0.05 and 15 cm for the Sun 
and Jupiter, respectively. The effect is much smaller for the Sun in spite of its much larger mass, due 
to its extremely slow motion in the barycentric frame. 

In view of the rapid decrease of gravitational deflection with increasing distance of closest ap- 
proach, it is extremely unlikely that a situation would arise where a VLBI observation would be close 
enough to a gravitating body for this correction to be of importance. In order to guard against such an 
unlikely situation, MASTERFIT code does perform the transit-time correction, however. To obtain 
the positions of the gravitational objects, we employ an iterative procedure, using the positions and 
velocities of the objects at signal reception time. If R(t r ) is the position of the gravitational object at 
signal reception time, t r , then that object’s position, R(t a ), at the time, t a , of closest approach of the 
ray path to the object was: 

R(t a ) = R(t r ) - \[t r - t a ] (2.29) 



c 


(2.30) 


We do this correction iteratively, using the velocity, V(£ r ), as an approximation of the mean velocity, 
V. Because v/c « 10 -4 , an iterative solution: 


Rn (ta) — R(tf) 


XM 

c 


|Rn-l(£a)| 


(2.31) 


rapidly converges to the required accuracy. 

Both the gravitational potential effects and the curved wave front effects are calculated inde- 
pendently of each other since the gravitational effects are a small perturbation (< l".75 in angle for 
Sun-grazing rays, or w 9 microradians). 
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2.2 TIME INFORMATION 


Before continuing the description of the geometric model, a few words must be said about time-tag 
information and the time units which will appear as arguments below. The epoch timing information 
in the data is taken hom the UTC time tags in the data stream at station #1. This time is converted 
to Terrestrial Dynamic Time (TDT) and is also used as an argument to obtain an a priori estimate 
of Earth orientation. The conversion consists of the following components: 

TDT = ( TDT - TAI) + (TAI - UTC BIH ) + ( UTCbih ~ UTC 0 ) 

+ (UTC 0 - UTC i ) + UTC i (2.32) 

where in seconds: 

TDT - TAI = 32.184 (2.33) 

and where TAI is atomic time. TAI - UTCbih = published integer second offset after 0 h , January 
1, 1972, and 

TAI - UTCbih = 9.8922417 + 3.0 x 10" 8 x ( UTC BIH - UTC 0 bih) (2.34) 

between 0 h , January 1, 1968, and 0 h , January 1, 1972. UTCbih — UTCo bih = number of UTC 
seconds relative to January 1, 1972. This is a negative number prior to that date. The software 
will not allow this quantity to be obtained prior to 1968. UTCbih ~ UTCo = the offset in UTC 
seconds between BIH UTC and the UTC clock at some secondary standard (usually NBS in Boulder 
for DSN observations). This can be obtained from BIH Circular D (typical reference is BIH Annual 
Report, 1981). In practice as of January, 1972, all that we do is use a linear interpolation between 
( UTCbih — UTChbs) data points as published in Circular D. The approximation usually is made 
that the clock at station #1 is very close to the NBS clock, e.g., UTCo — UTC\ < 5-10 ps. Since 
this time is used as epoch time in the observations, the major consequence resulting from an error in 
this assumption is to make an error in the estimation of UTl— UTC of one second per second of error 
in ( UTCi — UTCj). If short baselines are involved, then an error in epoch time causes an error of 
« Bu>e At = 7.3 x 10 -6 cm per km baseline per ps of clock error. This equals 0.0073 cm of baseline 
orientation error for a 10-^s clock epoch error on a 100-km baseline. 

A priori UTl— UTC and pole positions are normally obtained by interpolation of the BIH Circular 
D smoothed values. However, any other source of UTl— UTC and pole position could be used provided 
it is a function of UTC, and is expressed in a left-handed coordinate system (see section 2.6 below). 
Part of the documentation for any particular set of results should clearly state what were the values 
of UTl— UTC and pole position used in the data reduction process. 

For the Earth model based on the new IAU conventions, the following definitions are employed 
throughout (Kaplan, 1981): 

1. Julian date at epoch J2000 = 2451545.0. 

2. All time arguments denoted by T below are measured in Julian centuries of 36525 days of the 
appropriate time relative to the epoch J2000, i.e., T = (J D — 2451545. 0)/36525. 

3. For the time arguments used to obtain precession, nutation, or to reference the ephemeris, 
Barycentric Dynamic Time ( TDB ) is used. This is related to Terrestrial Dynamic Time ( TDT ) 
by the following: 


where 


TDB = TDT + 0. *001658 sin(<7 + 0.0167 sin(g)) 

(357.°528 + 35999. °050 TDT) x 2w 
9 ~ 360° 


(2.35) 

(2.36) 
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2.5 STATION LOCATION TIME DEPENDENCE 


The pre-1984 model considered the three coordinates of station i: r, Pi , A,-, Zi (radius off spin 
axis, longitude, and height above the equator, respectively) to be time- invariant. In investigations of 
tectonic motion, however, a new set of coordinates is usually solved for in the least-squares estimation 
process for each VLBI session. Post-processing software then makes linear fits to these results to infer 
bounds on the magnitude of the time rate of change of the station location. Care must be taken that 
the correlations of coordinates estimated at different epochs are accounted for properly. The advantage 
of this approach is that the contribution of each session to the overall slope may be independently 
evaluated, since it is clearly isolated. This procedure is somewhat inconvenient in practice, and an 
alternative is to introduce the time rates of change of the station coordinates as new parameters in 
MASTERFIT. The model is linear, with the coordinates at time t expressed as 


•Pi = r *p, + ’’•Pi (* — M 

(2.37) 

= A° + A,(t — t 0 ) 

(2.38) 

Zi = Zi +Zi(t- to) 

(2.39) 


Here to is a reference epoch, at which the station coordinates are (r° pj , A°, z°). 

As an alternative to linear time dependence of the station cylindrical coordinates, the MERIT 
standard model of tectonic plate rotation is optionally available in MASTERFIT. It is described in 
an addition to the MERIT standards document (Melbourne et al., 1985), and was denoted AMO-2 in 
the original paper (Minster and Jordan, 1978). Time dependence of the Cartesian station coordinates 
is expressed as 

*,• = x? + (wjz? - wjy?)(t - t 0 ) (2.40) 

Vi = Vi + [vUi ~ ~ to) (2.41) 

Zi = z° + (w 3 x y° - w£x°)(t - t 0 ) (2.42) 

where wj , are velocities of the plate j on which station i resides. Table I gives a list of the rotation 
rates for the 11 plates in the AMO-2 model. 


Table I 

Plate Rotation Velocities'): Minster- Jordan AMO-2 Model 


Plate 

w x 

Wy 

<*>* 

PCFC 

-0.12276 

0.31163 

-0.65537 


-0.63726 

-1.33142 

0.72556 

NAZC 

-0.09086 

-0.53281 

0.63061 

CARB 

-0.02787 

-0.05661 

0.10780 

SOAM 

-0.05604 

-0.10672 

-0.08642 

ANTA 

-0.05286 

-0.09492 

0.21570 

INDI 

0.48372 

0.25011 

0.43132 

AFRC 

0.05660 

-0.19249 

0.24016 

ARAB 

0.27885 

-0.16744 

0.37359 

EURA 

-0.03071 

-0.15865 

0.19605 

NOAM 

0.03299 

-0.22828 

-0.01427 


t units are /ideg/year 


At present, there is no facility in MASTERFIT to compute partial derivatives with respect to 
the plate velocities, or to solve for these quantities. 


i 
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2.4 TIDAL EFFECTS 


As an initial step in calculating the geometric delay, we need to consider the effects of crustal 
motions on station locations. Among these deformations are solid Earth tides, tectonic motions, and 
alterations of the Earth’s surface due to local geological, hydrological, and atmospheric processes. One 
possibility is to not model crustal movement other than that due to solid Earth tides, and allow the 
other effects to manifest themselves as temporal changes of the Earth-fixed baseline. Such a strategy 
corrupts the estimation of global orientation parameters from a finite set of baselines, and is a known 
weakness («1— 10 cm/year) of the simplified form of the current model. 

In the standard terrestrial coordinate system, tidal effects modify the station location ro by an 
amount 

A = Ajoj + A ocn ■+■ A a < m + A po j (2.43) 

where the four terms are due to solid Earth tides, ocean loading, atmosphere loading, and pole tide, 
respectively. Other Earth-fixed effects would be incorporated by augmenting the definition of A. 

2.4.1 Solid Earth Tides 


Alteration of the positions of the stations by solid Earth tides is rather complicated due to their 
coupling with the ocean tides, and the effects of local geology. We have chosen to gloss over these 
complications initially, and to incorporate the simple quadrupole response model described by J. G. 
Williams (1970), who used Melchior (1966) as a reference. Let r t be the Cartesian coordinates of a 
station in a geocentric frame (x axis vertical, y axis eastward, and the z axis northward) on a spherical 
Earth. F'urther, let R, be the position of a perturbing source in the terrestrial reference system, and 
To the station position in the same coordinates. If h and l are the Love numbers, rp a phase shift 
of the tidal effects, and r p the phase-shifted station location in the terrestrial system, then the tidal 
displacements are 


£ = ^2 [fcfih*, 1?2»» !?3.] T 


(2.44) 
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where 


9i. = 




[r„ ■ It,] 2 rlR\ 


92, = 


R b . 

3 ^* r P f- . » lfv - _ Y .. 1 l r pl 

R* 


[r p • R,](y.x p - X.y p \ 


\J x P + y\ 


93 • = 


3/t.r p 

-rT [1p 


r p • R, ] 


\J x l + yl z » — , Zp \ x p x , + y P Y . ] 


\] x l + y\ 


(2.45) 

(2.46) 


(2.47) 


H, is the ratio of the mass of the disturbing object, a, to the mass of the Earth, and 

R. = [X„Y„Z.] T (2.48) 

is the vector from the center of the Earth to that body. The summation is over all tide-producing 
bodies, of which we include only the Sun and the Moon. If the tidal effect at time ti is desired, the 
position of the tide-producing mass at time 

tx - St = ti - |R(ti -6t)\/c (2.49) 

should be used (a nuance we have not yet incorporated). 

The phase-shifted station vector above is calculated employing a phase lag, or equivalently, em- 
ploying a right-handed rotation, L, through an angle xf> about the z axis of date, r p = Lro- This lag 
matrix, L , is: 

( cos rp sin rp 0 \ 

— sin^ cos^ 0 I 

0 0 1 / 


12 


(2.50) 


By a positive value of if> we mean that the peak response on an Earth meridian occurs at a time 
St = r/>/uE after that meridian containing ro crosses the tide-producing object, where ug is the 
angular rotation rate of the Earth. In the vertical component, the peak response occurs when the 
meridian containing r p also includes R, . The difference between geodetic and geocentric latitude can 
affect this model on the order of the (tidal effect) /(flattening factor) as 0.1 cm. 

To convert this locally referenced strain, 6, to the Earth-fixed frame, two additional rotations 
must be performed. The first, W, rotates by an angle, — <f>,, about the y axis to an equatorial system. 
The second, V, rotates about the resultant s axis by angle, —A,, to bring the displacements into the 
standard geocentric coordinate system. The result is 


where 


and 


A aoi = VW6 


W = 


V = 



Actually, the product of these two matrices is coded: 


VW 


( cos A, cos 4>, — sin A, 
sin A, cos 4> t cos A, 
sin^, 0 


— cos A, sin <f> 

— sin A, sin <f> 

cos <f>. 



(2.51) 


(2.52) 


(2.53) 


(2.54) 


2.4.2 Ocean Loading 

This section is concerned with secondary tidal effects, the elastic response of the Earth’s 
crust to ocean tides, which move the observing stations to the extent of a few cm. Such effects are 
commonly labeled “ocean loading.” The most complete recent model appears to be that described by 
Pagiatakis (1982) and by Pagiatakis, Langley, and Vanicek (1982), which is implemented as an option 
in the current MASTERFIT. The formulation and coding are general enough, however, to permit 
other inputs to be used in place of the Pagiatakis ocean loading model. Because the station motions 
caused by response to ocean tides appear to be limited to approximately 3 cm for sites well removed 
from the coast, no estimation capability was deemed necessary at present. This decision is supported 
by the fact that for locations near the coast, where the effects may be more sizable, and which would 
thus be expected to produce data useful in parameter estimation, the elastic response modeling is as 
yet inadequate (Agnew, 1982). As suggested in section 8 of the initial version of this report (V1.0), 
local Earth motion can be partially accounted for by varying the Love numbers for each station. The 
present model entails deriving an expression for the locally referenced displacement 6 due to ocean 
loading. In the vertical, N-S, E-W local coordinate system at time t, 

N 

6j = ^2 g cos(w,t + Vi - Si) 

»=i 

The quantities w,- (frequency of tidal constituent t) and V, (astronomical argument of constituent t) 
depend only on the ephemeris information (positions of the Sun and Moon). On the other hand the 
amplitude £* and Greenwich phase lag Sj of each tidal component j are determined by the particular 
model assumed for the deformation of the Earth. As of December 1984, software for calculating these 
deformations at an arbitrary point on the Earth’s surface exists at JPL only for the Pagiatakis-Langley 
model. Six tidal components are included [IV = 6 in Eq. (2.55)]: the M 2 , S 2 , Ki, 0\, N 2 , and Pi tides, 


(2.55) 
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all of which have periods close to either 12 or 24 hours. The local displacement vector is transformed 
via Eqs. (2.54) and (2.51) to the displacement Aocn in the standard geocentric frame. 

Input to MASTERFIT provides for specification of an arbitrary number of frequencies and astro- 
nomical arguments and Vi, followed by tables of the local distortions and their phases, £? and S- t 
calculated from the ocean tidal loading model of choice. In particular, longer-period tidal constituents 
can be accommodated in this fashion. 

There are presently four choices of models. As mentioned above, the six components of the 
Pagiatakis-Langley model can be calculated by separate software for an arbitrarily located station 
to generate input tables for MASTERFIT. This has been done for all stations commonly employed 
in JPL VLBI experiments. There is no comparable facility to obtain amplitudes and phases for 
the Agnew (1982), Scherneck (1983), or Goad (1983) models. Consequently, for these three models, 
MASTERFIT input tables only exist for the limited set of stations considered by these authors. 
Agnew considers only five components (omitting Pi), while Scherneck includes five components in 
addition to the Pagiatakis-Langley set: K?, Q lt Mj , M m , and S. a - These all have amplitudes of 1 
mm or smaller and, thus, are not expected to be significant at the present level of VLBI experimental 
accuracy. Goad’s MERIT standard model only specifies vertical displacements, but includes three 
components (Jf 2> Qi and Mj) in addition to Pagiatakis-Langley. Due to their bulk, none of the tables 
of tidal amplitudes is reproduced here, but are available on request in computer-readable form. The 
default tidal model in MASTERFIT remains the Williams quadrupole solid Earth tide model with no 
ocean loading. 

2.4.8 Atmosphere Loading 

By analogy with the consequences of ocean tides that were considered in the previous section, 
a time-varying atmospheric pressure distribution can induce crustal deformation. A recent paper by 
Rabbel and Schuh (1986) estimates the effects of atmospheric loading on VLBI baseline determina- 
tions, and concludes that they may amount to many millimeters of seasonal variation. In contrast to 
ocean tidal effects, analysis of the situation in the atmospheric case does not benefit from the presence 
of a well-understood periodic driving force. Otherwise, estimation of atmospheric loading via Green’s 
function techniques is analogous to methods used to calculate ocean loading effects. Rabbel and Schuh 
recommend a simplified form of the dependence of the vertical crust displacement on pressure distri- 
bution. It involves only the instantaneous pressure at the site in question, and an average pressure 
over a circular region C of radius R = 2000 km surrounding the site. The expression for the vertical 
displacement (mm) is: 

A r = — 0.35p - 0.55p (2.56) 

where p is the local pressure anomaly, and p the pressure anomaly within the 2000-km circular region 
mentioned above (both quantities are in mbar). Note that the reference point for this displacement is 
the site location at standard pressure (1013 mbar). The locally referenced Ar is transformed to the 
standard geocentric coordinate system via the transformation (2.54). 

It was decided to incorporate this rudimentary model into MASTERFIT as an optional part of 
the model, with an additional mechanism for characterizing p. The two-dimensional surface pressure 
distribution surrounding a site is described by 

p(x, y) = Aq + A\x + A?y + A31 2 + A^xy + Agy 2 (2-57) 


where x and y are the local East and North distances of the point in question from the VLBI site. 
The pressure anomaly p may then be evaluated by the simple integration 



p = dx d y p( x > v) / f f c dx d y 

(2.58) 

giving 


p = Ao + (A3 + As)R /4 

(2.59) 


14 



It remains the task of the data analyst to perform a quadratic fit to the available weather data 
to determine the coefficients Aq-s- Future advances in understanding the atmosphere-crust elastic 
interaction can probably be accommodated by adjusting the coefficients in Eq. (2.56). 


2.4.4 Pole Tide 

In addition to the effects of ocean tides and atmospheric pressure variations, another secondary 
tidal effect is the displacement of a station by the elastic response of the Earth’s crust to shifts in the 
spin axis orientation. The spin axis is known to describe a circle of » 20-m diameter at the north pole. 
Depending on where the spin axis pierces the crust at the instant of a VLBI measurement, the “pole 
tide” displacement will vary from time to time. This effect must be included if centimeter accuracy 
is desired. 

Yoder (1984) derived an expression for the displacement of a point at latitude <f>, longitude A due 
to the pole tide: 


6 = — U>E ^ [sin <j> cos <p{x cos A + y sin A)hr 
+ cos 2 <j>(x cos A + y sin A)Z^ 

+ sin$(— xsinA + ycosA)/A] (2.60) 


Here we is the rotation rate of the Earth, R the radius of the (spherical) Earth, g the acceleration 
due to gravity at the Earth’s surface, and h and l the customary Love numbers. Displacements of the 
spin axis from the 1903.0 CIO pole position along the x and y axes are given by x and y. Eq. (2.60) 
shows how these map into station displacements along the unit vectors in the radial (r), latitude (^), 
and longitude (A) directions. With the standard values we = 7.292 x 10“ 6 rad/sec, R = 6378 km, 
and g = 980.665 g/cm 2 , the factor w\Rjg = 3.459 X 10 -3 . Since the maximum values of i and y are 
of the order of 10 meters, and h <=s 0.6, l « 0.08, the maximum displacement due to the pole tide is 1 
to 2 cm, depending on the location of the station (<^, A). 

As in the case of the previously considered tidal effects, the locally referenced displacement S is 
transformed via the transformation (2.54) to give the displacement A po i in the standard geocentric 
coordinate system. The pole tide effect has been coded as an optional part of the MASTERFIT model. 
It is only applied if specifically requested, t.e., the default model contains no pole tide contributions 
to the station locations. 

After each of the locally referenced tidal displacements has been transformed to standard terres- 
trial coordinates, the station location is 


Yt To A §0 { Aocn ZA a im A poi 


(2.61) 
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2.5 TRANSFORMATION FROM TERRESTRIAL TO CELESTIAL 
COORDINATE SYSTEMS 

The Earth is approximately an oblate spheroid, spinning in the presence of two massive moving 
objects (the Sun and the Moon) which are positioned such that their time-varying gravitational effects 
not only produce tides on the Earth, but also subject it to torques. In addition, the Earth is covered 
by a complicated fluid layer, and also is not perfectly solid internally. As a result, the orientation 
of the Earth is a very complicated function of time, which to first order can be represented as the 
composite of a time-varying rotation rate, a wobble, a nutation, and a precession. The exchange of 
angular momentum between the solid Earth and the fluids on its surface is not readily predictable, 
and thus must be continually determined experimentally. Nutation and precession are well modeled 
theoretically. However, at the accuracy with which VLBI can determine baseline vectors, even these 
models are not completely adequate. 

Currently, the rotational transformation, Q, of coordinate frames from the terrestrial frame to 
the celestial geocentric frame is composed of 6 separate rotations (actually 10, since the nutation 
transformation, N, consists of 3 transformations in itself, as does the “perturbation” transformation, 
A) applied to a vector in the terrestrial system: 

Q = CIPNUXY (2.62) 

Thus, if r t is a station location expressed in the terrestrial system, e.g., the result of (2.61), that 
location, r c , expressed in the celestial system is 


r c = Qr t (2.63) 

This particular formulation follows the historical path of astrometry, and is couched in that lan- 
guage. While aesthetically unsatisfactory with modern measurement techniques, such a formulation is 
currently practical for intercomparison of techniques and for effecting a smooth inclusion of the inter- 
ferometer data into the long historical record of astrometric data. Much more pleasing aesthetically 
would be the separation of Q into two rotation matrices: 

Q = Q1Q2 (2-64) 

where Q2 are those rotations to which the Earth would be subjected if all external torques were 
removed (approximately UXY above), and where Qi are those rotations arising from external torques 
(approximately CIPN above). Even then, the tidal response of the Earth prevents such a separation 
from being perfectly realized. Eventually, the entire problem of obtaining the matrix Q, and the tidal 
effects on station locations should be done numerically. Note that since we rotate the Earth rather 
than the celestial sphere, our rotation matrices, O, P, and N, will be the transposes of those used to 
rotate the celestial system of J2000 to a celestial system of date. 

2.6 UT1 AND POLAR MOTION 

The first transformation, Y, is a right-handed rotation about the x axis of the terrestrial frame 
by an angle © 2 - Currently, the terrestrial frame is the 1903.0 CIO frame, except that the positive y 
axis is at 90 degrees east (Moscow). The x axis is coincident with the 1903.0 meridian of Greenwich, 
and the 1 axis is the 1903.0 standard pole. 

/! 0 .0 \ 

Y — I 0 cos ©2 sin ©2 I (2.65) 

1^0 —sin ©2 cos ©2 J 

where ©2 is identical to the B y that may be readily obtained from the pole position published by 

BIH. 
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The next rotation in sequence is the right-hand rotation through an angle 0i about the y axis 
obtained after the previous rotation has been applied: 

( cos 0i 0 —sin 0i \ 

0 10 
sin 0i 0 cos 0i J 

In this rotation, 0i is identical to the 0 Z obtainable from the BIH value for the pole position. Note 
that we have incorporated in the matrix definitions the transformation from the left-handed system 
used by BIH to the right-handed system we use. Note also that instead of BIH data used as a pole 
definition, we could instead use any other source of polar motion data provided it was represented in 
a left-handed system. The only effect would be a change in the definition of the terrestrial reference 
system. 

The application of U XY ” to a vector in the terrestrial system of coordinates expresses that vector 
as it would be observed in a coordinate frame whose z axis was along the Earth’s ephemeris pole. 
The third rotation, U, is about the resultant z axis obtained by applying “ XY ”. It is a rotation 
through the angle, — H, where H is the hour angle of the true equinox of date (t.e., the dihedral angle 
measured westward between the xz plane defined above and the meridian plane containing the true 
equinox of date). The equinox of date is the point defined on the celestial equator by the intersection 
of the mean ecliptic with that equator. It is that intersection where the mean ecliptic rises from below 
the equator to above it (ascending node). 


( 2 . 66 ) 


( cos H — sin H 0 
sin H cos H 0 

0 0 1 

This angle H is composed of two parts: 

H = h y + as 


(2.67) 


( 2 . 68 ) 


where h n is the hour angle of the mean equinox of date, and as is the difference in hour angle of 
the true equinox of date and the mean equinox of date, a difference which is due to the nutation of 
the Earth. This set of definitions is cumbersome and couples the nutation and precession effects into 
Earth rotation measurements. However, in order to provide a direct estimate of conventional UT 1 it 
is convenient to endure this historical approach, at least for the near future. 

UT 1 (universal time) is defined to be such that the hour angle of the mean equinox of date is 
given by the following expression (Aoki et al., 1982, and Kaplan, 1981): 

= UT1 + 6 h 41 m 50*. 54841 + 8640184*.812866 T„ 

+ 0* .093104 T% - 6*. 2 X 10~ 6 T* (2.69) 


where 


T u 


(Julian UTl date) - 2451545.0 


36525 


(2.70) 


The actual equivalent expression which is coded is: 


=2n(UTl Julian day fraction) + 67310*. 54841 

+ 8640184*. 812866 T u + 0* .093104 T? - 6*. 2 X 1(T 6 T® 


(2.71) 


This expression produces a time, UTl, which tracks the Greenwich hour angle of the real Sun to 
within 16 m . However, it really is sidereal time, modified to fit our intuitive desire to have the Sun 
directly overhead at noon on the Greenwich meridian. Historically, differences of UTl from a uniform 
measure of time, such as atomic time, have been used in specifying the orientation of the Earth. Note 
that this definition has buried in it the precession constant since it refers to the mean equinox of date. 
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By the very definition of "mean of date” and “true of date” , natation causes a difference in the 
hour angles of the mean equinox of date and the true equinox of date. This difference, called the 
“equation of equinoxes”, is denoted by as and is obtained accordingly: 


where the vector 



(2.72) 


(2.73) 


is the unit vector, in true equatorial coordinates of date, toward the mean equinox of date. In mean 
equatorial coordinates of date, this same unit vector is just (1,0, 0) T . The matrix N~ l is just the 
inverse (or equally, the transpose) of the transformation matrix N which will be defined below to 
effect the transformation from true equatorial coordinates of date to mean equatorial coordinates of 
date. 

Depending on the smoothing used to produce the a priori UTl — UTC series, the short-period 
(t < 35 days) fluctuations in UTl due to changes in the latitude and size of the mean tidal bulge may 
or may not be smoothed out. Since we want as accurate an a priori as possible, it may be necessary 
to add this effect to the UTl a priori obtained from the series, UTl, moothed . If this option is selected, 
then the desired a priori UTl is given by 


UTla priori — UTl, m oothed "b AZ/TT 


(2.74) 


UTl, moo thed represents an appropriately smoothed o priori measurement of the orientation of the 
Earth (t'.e., typically BIH Circular D smoothed or, even better, UT1R), for which the short period 
(t < 35 days) tidal effects have either been averaged to zero, or, as in the case of UT1R, removed 
before smoothing. This AZ7T1 can be represented as 


AUT1 = 



(2.75) 


where N is chosen to include all terms with a period less than 35 days. There are no other contri- 
butions until a period of 90 days is reached. However, these long-period terms are included by the 
measurements of the current Earth-orientation measurement services. The values for ki 3 and Ai, along 
with the period involved, are given in Table II. The a< for t = 1, 5 are just the angles defined below 
in the nutation series as l, l', F, D, and fl, respectively. In Table II, the sign of the 14.73 day term 
has been changed from that published [Yoder (1982)] to remove a sign error in the article referenced. 
Otherwise, the values in Table II are those given by BIH (1982) as derived from Yoder et al. (1981). 

It might be appropriate at this point to describe the interpolation method used in MASTERFIT 
to obtain a priori polar motion and UTl values. These are normally available as tables at 5-day 
intervals, from either BIH or the IRIS project (IAG, 1986). Linear interpolation is performed for 
all three quantities. If the short-period tidal terms At/Tl are present in the tabular values, they 
are subtracted before interpolation, and added back to the final value. With the present accuracy of 
determinations of pole position and Z7T1, linear interpolation over a 5-day interval may be inadequate, 
possibly giving rise to 0.1 to 0.2 ms errors in UTl. Quadratic spline interpolation is being considered 
as an alternative. Even with the present code, however, the highest possible accuracy may be achieved 
by performing the interpolation externally to MASTERFIT, and supplying it with tables of values 
more closely spaced in time for the final internal linear interpolation. The Kalman-filtered UTPM 
values of Eubanks et al. (1984) are ideally suited for this purpose. 
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Table II 

Periodic Tidally Induced Variations in UT1 
with Periods Less than 35 Days 


Index 

i 

Period 

(days) 

Argument coefficient 
^»1 ^*2 ^*3 ^*4 ^*6 

A, 

(O'.OOOl) 

i 

5.64 

l 

0 

2 

2 

2 

-0.02 

2 

6.85 

2 

0 

2 

0 

1 

-0.04 

3 

6.86 

2 

0 

2 

0 

2 

-0.10 

4 

7.09 

0 

0 

2 

2 

1 

-0.05 

5 

7.10 

0 

0 

2 

2 

2 

-0.12 

6 

9.11 

1 

0 

2 

0 

0 

-0.04 

7 

9.12 

1 

0 

2 

0 

1 

-0.41 

8 

9.13 

1 

0 

2 

0 

2 

-0.99 

9 

9.18 

3 

0 

0 

0 

0 

-0.02 

10 

9.54 

-1 

0 

2 

2 

1 

-0.08 

11 

9.56 

-1 

0 

2 

2 

2 

-0.20 

12 

9.61 

1 

0 

0 

2 

0 

-0.08 

13 

12.81 

2 

0 

2 

-2 

2 

0.02 

14 

13.17 

0 

1 

2 

0 

2 

0.03 

15 

13.61 

0 

0 

2 

0 

0 

-0.30 

16 

13.63 

0 

0 

2 

0 

1 

-3.21 

17 

13.66 

0 

0 

2 

0 

2 

-7.76 

18 

13.75 

2 

0 

0 

0 

-1 

0.02 

19 

13.78 

2 

0 

0 

0 

0 

-0.34 

20 

13.81 

2 

0 

0 

0 

1 

0.02 

21 

14.19 

0 

-1 

2 

0 

2 

-0.02 

22 

14.73 

0 

0 

0 

2 

-1 

0.05 

23 

14.77 

0 

0 

0 

2 

0 

-0.73 

24 

14.80 

0 

0 

0 

2 

1 

-0.05 

25 

15.39 

0 

-1 

0 

2 

0 

-0.05 

26 

23.86 

1 

0 

2 

-2 

1 

0.05 

27 

23.94 

1 

0 

2 

-2 

2 

0.10 

28 

25.62 

1 

1 

0 

0 

0 

0.04 

29 

26.88 

■1 

0 

2 

0 

0 

0.05 

30 

26.98 

-1 

0 

2 

0 

1 

0.18 

31 

27.09 

-1 

0 

2 

0 

2 

0.44 

32 

27.44 

1 

0 

0 

0 

-1 

0.53 

33 

27.56 

1 

0 

0 

0 

0 

-8.26 

34 

27.67 

1 

0 

0 

0 

1 

0.54 

35 

29.53 

0 

0 

0 

1 

0 

0.05 

36 

29.80 

1 

-1 

0 

0 

0 

-0.06 

37 

31.66 

-1 

0 

0 

2 

-1 

0.12 

38 

31.81 

-1 

0 

0 

2 

0 

-1.82 

39 

31.96 

-1 

0 

0 

2 

1 

0.13 

40 

32.61 

1 

0 

-2 

2 

-1 

0.02 

41 

34.85 

-1 

-1 

0 

2 

0 

-0.09 
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It is convenient to apply “UXY ” as a group. To parts in 10 12 , XY = YX. However, with the 
same accuracy UXY ± XYU . Neglecting terms of 0 (0 2 ) (which produce station location errors of 
approximately 6 X 10 -6 meters): 


( cos H 
sin H 
sin 0i 


— sin H 
cos H 

— sin 02 


sin ©i cos. 


(2.76) 


2.7 NUTATION 

With the completion of the UTl and polar motion transformations, we are left with a station 
location vector, r date- This is the station location relative to true equatorial celestial coordinates of 
date. The last set of transformations are nutation, N, precession, P, and the perturbation rotation, 
Cl, applied in that order. These transformations give the station location, r c , in celestial equatorial 
coordinates: 

r c = CIPNTdau (2.77) 

The transformation matrix N is a composite of three separate rotations (Melbourne et al., 1968): 

1. A(e): true equatorial coordinates of date to ecliptic coordinates of date. 

(10 0 \ 

A(e) = I 0 cose sine | (2.78) 

y 0 — sin e cos e J 

2. (^(Sip): nutation in longitude from ecliptic coordinates of date to mean ecliptic coordinates of 
date. 


C t (6iP) = 


( cos Sip 
— sin Sip 
0 


sin Sip 
cos Sip 
0 



where Sip is the nutation in ecliptic longitude. 


(2.79) 


3. A T (e): ecliptic coordinates of date to mean equatorial coordinates. 

In ecliptic coordinates of date, the mean equinox is at an angle Sip = tan -1 (y^-/a:y). Se = e — e 
is the nutation in obliquity, and e is the mean obliquity (the dihedral angle between the plane of the 
ecliptic and the mean plane of the equator). “Mean” as used in this section implies that the short- 
period (T < 18.6 years) effects of nutation have been removed. Actually, the separation between 
nutation and precession is rather arbitrary, but historical. The composite rotation is: 

N = A T (e)C T {Sip)A(e) 

( cos Sip cos s sin Sip sins sin Sip 

— cosssi nSip cosscosecos Sip + sin s sin e cos e sin e cos Sip — sin e cos e 

— sine sin Srp sine cose cob Sip — cose sine sine sin ecos Sip + cose cose 

The 1980 IAU nutation model (Seidelmann, 1982, and Kaplan, 1981) is used to obtain .the values 
for Sip and e — e. The mean obliquity is obtained from Lieske et al. (1977) or from Kaplan (1981): 

e = 23° 26' 21."448 - 46."8150 T - 5."9 X 10~ 4 T 2 + l."813 X lO -3 ! 0 (2.81) 

This nutation in longitude (Sip) and in obliquity ( Se = e — e ) can be represented by a series expansion 

of the sines and cosines of linear combinations of five fundamental arguments. These are (Kaplan, 
1981, Cannon, 1981): 


(2.80) 
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1. the mean anomaly of the Moon: 


ai = l= 485866".733 + (1325 r + 715922".633) T 

+ 31".310T 2 + 0".064 T 3 (2.82) 

2. the mean anomaly of the Sun: 

a 2 = 1'= 1287099". 804 + (99 r + 129258l".224) T 

- 0".577 T 2 - 0."012 T 3 (2.83) 

3. the mean argument of latitude of the Moon: 

a 3 = F= 335778".877 + (1342 r + 295263".137) T 

- 13". 257 T 2 + 0".011 T 3 (2.84) 

4. the mean elongation of the Moon from the Sun: 

a 4 = D = 1072261".307 + (1236 r + 110560l".328) T 

- 6". 891 T 2 + 0".019 T 3 (2.85) 

5. the mean longitude of the ascending lunar node: 

a 5 = A = 450160".280 - (5 r + 482890".539) T 

+ 7". 455 T 2 + 0".008 T 3 (2.86) 

where l r = 360° = 1296000". 

With these fundamental arguments, the nutation quantities then can be represented by 


and 


N 


i=i L 


(v4oy + AijT) sin 


N 


1=1 L 


(Boj + BijT) cos 



(2.87) 


( 2 . 88 ) 


where the various values of a<, kji, Ay, and Bj are tabulated in Table III. 

An additional set of terms can be optionally added to the nutations Sip and Se in Eqs. (2.87) and 
(2.88). These include the out-of-phase nutations, and the free-core nutations (Yoder, 1983) with period 
Uf (nominally 460 days). The “nutation tweaks* A if> and As are arbitrary constant increments of 
the nutation angles Sip and Se. Unlike the usual nutation expressions, they have no time dependence. 
The out-of-phase nutations, which are not included in the IAU 1980 nutation series, are identical to 
Eqs. (2.87) and (2.88), with the replacements sin *-* cos: 


and 


N 

• 

■ 5 ' 

sip = ^2 

(Aay + A3 3 ‘T) cos 


1=1 

. 

.*=i j. 


AT 




j=il 


(B 2 j + B 3j T) sin 




U=i 


(2.89) 


(2.90) 
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Table III 


1980 IAU Theory of Nutation 


Index 

j 

Period 

( days ) 

k H 

Argument coefficient 
kj 2 fcys fc/4 k 

>6 

Aoy Aiy 

( 0 ". 0001 ) 

Bo ,- Bij 

( 0 ". 0001 ) 

1 

6798.4 

0 

0 

0 

0 

1 

-171996 

- 174.2 

92025 

8.9 

2 

3399.2 

0 

0 

0 

0 

2 

2062 

0.2 

-895 

0.5 

3 

1305.5 

-2 

0 

2 

0 

1 

46 

0.0 

-24 

0.0 

4 

1095.2 

2 

0 

-2 

0 

0 

11 

0.0 

0 

0.0 

5 

1615.7 

-2 

0 

2 

0 

2 

-3 

0.0 

1 

0.0 

6 

3232.9 

1 

-1 

0 

-1 

0 

-3 

0.0 

0 

0.0 

7 

6786.3 

0 

-2 

2 

-2 

1 

-2 

0.0 

1 

0.0 

8 

943.2 

2 

0 

-2 

0 

1 

1 

0.0 

0 

0.0 

9 

182.6 

0 

0 

2 

-2 

2 

-13187 

- 1.6 

5736 

- 3.1 

10 

365.3 

0 

1 

0 

0 

0 

1426 

- 3.4 

54 

- 0.1 

11 

121.7 

0 

1 

2 

-2 

2 

-517 

1.2 

224 

- 0.6 

12 

365.2 

0 

-1 

2 

-2 

2 

217 

- 0.5 

-95 

0.3 

13 

177.8 

0 

0 

2 

-2 

1 

129 

0.1 

-70 

0.0 

14 

205.9 

2 

0 

0 

-2 

0 

48 

0.0 

1 

0.0 

15 

173.3 

0 

0 

2 

-2 

0 

-22 

0.0 

0 

0.0 

16 

182.6 

0 

2 

0 

0 

0 

17 

- 0.1 

0 

0.0 

17 

386.0 

0 

1 

0 

0 

1 

-15 

0.0 

9 

0.0 

18 

91.3 

0 

2 

2 

-2 

2 

-16 

0.1 

7 

0.0 

19 

346.6 

0 

-1 

0 


1 

-12 

0.0 

6 

0.0 

20 

199.8 

-2 

0 

0 

2 

1 

-6 

0.0 

3 

0.0 

21 

346.6 

0 

-1 

2 

-2 

1 

-5 

0.0 

3 

0.0 

22 

212.3 

2 

0 

0 

-2 

1 

4 

0.0 

-2 

0.0 

23 

119.6 

0 

1 

2 

-2 

1 

4 

0.0 

-2 

0.0 

24 

411.8 

1 


0 

-1 

0 

-4 

0.0 

0 

0.0 

25 

131.7 

2 

1 


-2 


1 

0.0 

0 

0.0 

26 

169.0 



-2 

2 

1 

1 

0.0 

0 

0.0 

27 

329.8 


1 

-2 

2 

0 

-1 

0.0 

0 

0.0 

28 

409.2 

0 

1 

0 

0 

2 

1 

0.0 

0 

0.0 

29 

388.3 

-1 

0 

0 

1 

1 

1 

0.0 

0 

0.0 

30 

117.5 

0 

1 

2 

-2 

0 

-1 

0.0 

0 

0.0 

31 

13.7 

0 

0 

2 

0 

2 

-2274 

- 0.2 

977 

- 0.5 

32 

27.6 

1 

0 

0 

0 

0 

712 

0.1 

-7 

0.0 

33 

13.6 

0 

0 

2 

0 

1 

-386 

- 0.4 

200 

0.0 

34 

9.1 

1 

0 

2 

0 

2 

-301 

0.0 

129 

- 0.1 

35 

31.8 

1 

0 

0 

-2 

0 

-158 

0.0 

-1 

0.0 

36 

27.1 

-1 

0 

2 

0 

2 

123 

0.0 

-53 

0.0 

37 

14.8 

0 

0 

0 

2 

0 

63 

0.0 

-2 

0.0 

38 

27.7 

1 

0 

0 

0 

1 

63 

0.1 

-33 

0.0 

39 

27.4 

-1 

0 

0 

0 

1 

-58 

- 0.1 

32 

0.0 

40 

9.6 

-1 

0 

2 

2 

2 

-59 

0.0 

26 

0.0 
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Table III cont. 


1980 IAU Theory of Nutation 


Index 

j 

Period 

(days) 

kji 

Argument coefficient 
k-)2 k)3 

fcys 

Aoy Aiy 

(0".0001) 

Boy Bij 

(0".0001) 

41 

9.1 

i 

0 

2 

0 

1 

-51 

0.0 

27 

0.0 

42 

7.1 

0 

0 

2 

2 

2 

-38 

0.0 

16 

0.0 

43 

13.8 

2 

0 

0 

0 

0 

29 

0.0 

-1 

0.0 

44 

23.9 

1 

0 

2 

-2 

2 

29 

0.0 

-12 

0.0 

45 

6.9 

2 

0 

2 

0 

2 

-31 

0.0 

13 

0.0 

46 

13.6 

0 

0 

2 

0 

0 

26 

0.0 

-1 

0.0 

47 

27.0 

-1 

0 

2 

0 

1 

21 

0.0 

-10 

0.0 

48 

32.0 

-1 

0 

0 

2 

mm 

16 

0.0 

-8 

0.0 

49 

31.7 

1 

0 

0 

-2 

n 

-13 

0.0 

7 

0.0 

50 

9.5 

-1 

0 

2 

2 

■9 

-10 

0.0 

5 

0.0 

51 

34.8 

1 

1 

0 

-2 

0 

-7 

0.0 

0 

0.0 

52 

13.2 

0 

1 

2 

0 

2 

7 

0.0 

-3 

0.0 

53 

14.2 

0 

-1 

2 

0 

2 

-7 

0.0 

3 

0.0 

54 

5.6 

1 

0 

2 

2 

2 

-8 

0.0 

3 

0.0 

55 

9.6 

1 

0 

0 

2 

0 

6 

0.0 

0 

0.0 

56 

12.8 

2 

0 

2 

-2 

2 

6 

0.0 

-3 

0.0 

57 

14.8 

0 

0 

0 

2 

1 

-6 

0.0 

3 

0.0 

58 

7.1 

0 

0 

2 

2 

1 

-7 

0.0 

3 

0.0 

59 

23.9 

1 

0 

2 

-2 

1 

6 

0.0 

-3 

0.0 

60 

14.7 

0 

0 

0 

-2 

1 

-5 

0.0 

3 

0.0 

61 

29.8 

1 

-1 

0 

0 

0 

5 

0.0 

0 

0.0 

62 

6.9 

2 

0 

2 

0 

1 

-5 

0.0 

3 

0.0 

63 

15.4 

0 

1 

0 

-2 

0 

-4 

0.0 

0 

0.0 

64 

26.9 

1 

0 

-2 

0 

0 

4 

0.0 

0 

0.0 

65 

29.5 

0 

0 

0 

1 

0 

-4 

0.0 

0 

0.0 

66 

25.6 

1 

1 

0 

0 

0 

-3 

0.0 

0 

0.0 

67 

9.1 

1 

0 

2 

0 

0 

3 

0.0 

0 

0.0 

68 

9.4 

1 

-1 

2 

0 

2 

-3 

0.0 

1 

0.0 

69 

9.8 

-1 

-1 

2 

2 

2 

-3 

0.0 

1 

0.0 

70 

13.7 

-2 

0 

0 

0 

1 

-2 

0.0 

1 

0.0 

71 

5.5 

3 

0 

2 

0 

2 

-3 

0.0 

1 

0.0 

72 

7.2 

0 

-1 

2 

2 

2 

-3 

0.0 

1 

0.0 

73 

8.9 

1 

1 

2 

0 

2 

2 

0.0 

-1 

0.0 

74 

32.6 

-1 

0 

2 

-2 

1 

-2 

0.0 

1 

0.0 

75 

13.8 

2 

0 

0 

0 

1 

2 

0.0 

-1 

0.0 

76 

27.8 

1 

0 

0 

0 

2 

-2 

0.0 

1 

0.0 

77 

9.2 

3 

0 

0 

0 

0 

2 

0.0 

0 

0.0 

78 

9.3 

0 

0 

2 

1 

2 

2 

0.0 

-1 

0.0 

79 

27.3 

-1 

0 

0 

0 

2 

1 

0.0 

-1 

0.0 

80 

10.1 

1 

0 

0 

•4 

0 

-1 

0.0 

0 

0.0 











Table III cont. 


1980 IAU Theory of Nutation 


Index 

j 

Period 

(days) 

Argument coefficient 
kji kj2 kj3 

fcys 

Aoj A\j 

(0".0001) 

Boj Bij 

(0".0001) 

81 

14.6 

-2 

0 

2 

2 

2 

1 

0.0 

-1 

0.0 

82 

5.8 

-1 

0 

2 

4 

2 

-2 

0.0 

1 

0.0 

83 

15.9 

2 

0 

0 

-4 

0 

-1 

0.0 

0 

0.0 

84 

22.5 

1 

1 

2 

-2 

2 

1 

0.0 

-1 

0.0 

85 

5.6 

1 

0 

2 

2 

1 

-1 

0.0 

1 

0.0 

86 

7.3 

-2 

0 

2 

4 

2 

-1 

0.0 

1 

0.0 

87 

9.1 

-1 

0 

4 

0 

2 

1 

0.0 

0 

0.0 

88 

29.3 

1 

-1 

0 

-2 

0 

1 

0.0 

0 

0.0 

89 

12.8 

2 

0 

2 

-2 

1 

1 

0.0 

-1 

0.0 

90 

4.7 

2 

0 

2 

2 

2 

-1 

0.0 

0 

0.0 

91 

9.6 

1 

0 

0 

2 

1 

-1 

0.0 

0 

0.0 

92 

12.7 

0 

0 

4 

-2 

2 

1 

0.0 

0 

0.0 

93 

8.7 

3 

0 

2 

-2 

2 

1 

0.0 

0 

0.0 

94 

23.8 

1 

0 

2 

-2 

0 

-1 

0.0 

0 

0.0 

95 

13.1 

0 

1 

2 

0 

1 

1 

0.0 

0 

0.0 

96 

35.0 

-1 

-1 

0 

2 

1 

1 

0.0 

0 

0.0 

97 

13.6 

0 

0 

-2 

0 

1 

-1 

0.0 

0 

0.0 

98 

25.4 

0 

0 

2 

-1 

2 

-1 

0.0 

0 

0.0 

99 

14.2 

0 

1 

0 

2 

0 

-1 

0.0 

0 

0.0 

100 

9.5 

1 

0 

-2 

-2 

0 

-1 

0.0 

0 

0.0 

101 

14.2 

0 

-1 

2 

0 

1 

-1 

0.0 

0 

0.0 

102 

34.7 

1 

1 

0 

-2 

1 

-1 

0.0 

0 

0.0 

103 

32.8 

1 

0 

-2 

2 

0 

-1 

0.0 

0 

0.0 

104 

7.1 

2 

0 

0 

2 

0 

1 

0.0 

0 

0.0 

105 

4.8 

0 

0 

2 

4 

2 

-1 

0.0 

0 

0.0 

106 

27.3 

0 

1 

0 

1 

0 

1 

0.0 

0 

0.0 


Expressions similar to these are adopted for the free-core nutations: 

Sip = (A 0 o + -dioT) sin(w/T) + (A20 + A30 T) cos (w/T) (2.91) 

and 

Se = (Boo + BiqT) cos (w/T) + (B 20 + Bo oT) sin(w/T) (2.92) 

The nutation model thus contains a total of 856 parameters: Ai } - (t=0,3; j'= 1,106) and Bij (*=0,3; 
y=l,106) plus the free-nutation amplitudes A<o (*=0,3), Bio (t=0,3). The only nonzero a priori 
amplitudes are the Aoy, A iy, Bqj, Bij (j= 1,106) given in Table III. 

The nutation tweaks are just constant additive factors to the angles Sip and Se: 

Sip — Sip + Aip (2.93) 

and 

Se — Se + As (2.94) 
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As a temporary measure to provide for a better nutation model than the one in Table III, which 
is known to be deficient at the milliarcsecond level, there is now a MASTERFIT option to employ 
the annual and semiannual amplitudes of Herring et al. (1986). In terms of the present notation, 
and in the units of Table III, these revised amplitudes are: Ao ,9 = -13172.2, 80,9 = 5732.8, Ao.io = 
1471.0, Ho, io = 72.1 for the in-phase nutations, and Aj.o = -8.3, # 2,9 = -2.9, A. 2,10 = 15.8, . 82,10 = 
-2.2 for the out-of-phase terms. It is emphasised that, for the present, the default nutation model in 
MASTERFIT is just the 1980 IAU nutation model given in Table III. 


2.8 PRECESSION 

The next transformation in going from the terrestrial frame to the celestial frame is the rotation 
P. This is the precession transformation from mean equatorial coordinates of date to the equatorial 
coordinates of the reference epoch (e.g., J2000). It is the transpose of the matrix given on page 7 of 
Melbourne et al. (1968): 


Pn = cos £4 cos 9 ^ cos Z A — sin £4 sin Za ( 2 . 95 ) 

Pi 2 = cos cos 0,1 sin Za + sin £4 cos Za ( 2 . 96 ) 

P13 = cos $A sin ©A ( 2 . 97 ) 

P 21 = — sin $a cos ©a cos Za — cos £4 sin Za ( 2 . 98 ) 

822 = — sin £a cos ©a sin Za + cos £4 cos Za ( 2 . 99 ) 

P23 = — sin fA sin ©a (2.100) 

P31 = —sin © a cos Za (2.101) 

P32 = — sin ©a sin Za (2.102) 

P33 = cos ©a ( 2 . 103 ) 

With the angular units in arc seconds, the arguments are: 

g A = 2306.2181 T + 0.30188 T 2 + 0.017998 T 3 (2.104) 

Za = 2306.2181 T + 1.09468 T 2 + 0.018203 T 3 (2.105) 

©A = 2004.3109 T - 0.42665 T 2 - 0.041833 T 3 (2.106) 


These expressions are given by Lieske et al. (1977) and by Kaplan (1981). This completes the standard 
model for the orientation of the Earth. 

2.9 PERTURBATION ROTATION 

This standard model for the rotation of the Earth as a whole may need a small incremental 
rotation about any one of the resulting axes. Define this perturbation rotation matrix as 

0 = A x A y A z (2.107) 

where 

/! 0 0 \ 

A x = 0 1 -5©x (2.108) 

Vo se x 1 J 

with SO x being a small angle rotation about the x axis, in the sense of carrying y into z; 

/ 1 0 se y \ 

A y = 0 10 (2.109) 

V-6© y 0 1 J 
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with 56 w being a small angle rotation about the y axis, in the sense of carrying s into x; and 


( 1 -56, 0\ 

A, = I 56, 1 0 (2.110) 

V o 0 1 J 

with 56, being a small angle rotation about the s axis, in the sense of carrying x into y. For angles 
of the order of 1 arc second we can neglect terms of order SQ^Re as they give effects on the order of 
0.015 cm. Thus, in that approximation 


( 1 -56, 56„ \ 

fl = 56, 1 —56, (2.1H) 

V-56 w 56, 1 J 

In general, 

56, = 56<(t) = 6e i0 + SQiT + fi{T) (2.112) 

which is the sum of an offset, a time-linear rate, and some higher order or oscillatory terms. Currently, 
only the offset and linear rate are implemented. In particular, a non-sero value of 56„ is equivalent 
to a change in the precession constant. Setting 

56, = 56y = 56, = 0 (2.113) 

gives the effect of applying only the standard rotation matrices. 

Starting with the Earth-fixed vector, ro, we have in sections 2.3 through 2.8 above shown how 
we obtain the same vector, r c , expressed in the celestial frame: 

t c = (IPNUXY(t 0 + A) (2.114) 


2.10 EARTH ORBITAL MOTION 

We now wish to transform these station locations from a geocentric celestial reference frame 
moving with the Earth to a celestial reference frame which is at rest relative to the center of mass 
of the Solar System. In this Solar System barycentric frame we will use these station locations to 
calculate the geometric delay (see Section 2.1 above). We will transform this time interval so obtained 
back to the frame in which the time delay is actually measured by the interferometer — the frame 
moving with the Earth. 

Let E' be a geocentric frame moving with vector velocity = fic relative to a frame, E, at rest 
relative to the Solar System center of mass. Further, let r(t) be the position of a point (e.g., station 
location) in space as a function of time, t, as measured in the E (Solar System barycentric) frame. In 
the E' (geocentric) frame, there is a corresponding position r'(t') as a function of time, t'. We normally 
observe and model r'(f') as shown in sections 2.3 through 2.8 above. However, in order to calculate 
the geometric delay in the Solar System barycentric frame (E), we will need the transformations of 
r(t) and r'(t'), as well as of t and t', as we shift frames of reference. Measuring positions in units of 
light travel time, we have from Jackson (1975): 

r'(t') = r(t) + (7 - l)r(t) (2.115) 

P 

t' = 7 [t-r(t) -fi] (2.116) 

and for the inverse transformation: 

r(t) = r '(?) + ( 7 - l)r'(t') ■ ^ + 7 fit' (2.117) 

P 
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where 


t = 7 [t'+r'(t')fl 


(2.118) 


7 = (1 ~fi 2 ) ^ (2-119) 

Let t\ represent the time measured in the Solar System barycentric frame (E), at which a wave 
front crosses antenna 1 at position ri(ti). Let r2(ti) be the position of antenna 2 at this same time 
as measured in the Solar System barycentric frame. Also, let t% be the time measured in this frame at 
which that same wave front intersects station 2. This occurs at the position ^(fj)- Following section 
2.1 above, we can calculate the geometric delay — tj. Transforming this time interval back to the 
S' (geocentric) frame, we obtain 

t;’ - ti = 7(t; - ti) - 7[r 9 (tS) - ri(tx)] • P (2.120) 


Assume further that the motion of station #2 is rectilinear over this time interval. That assumption 
is not strictly true, but as discussed below it is sufficiently good such that the error made as a result 
of that assumption is much less than 1 cm in calculated delay. Thus, 


*■3(^2) = r a(ii) + Pai^ — *x) 


( 2 . 121 ) 


which gives: 


r 2( 4 2) — r x(*x) — r 2(*i) “ r x(*x) + Pafta — *i) 


and 


* 2 ' “ *x = 7(*2 “ *x) - 7[r 2 (fi) - ri(*i)] P ~ iPa P[*a ~ *x] 
= Tf [(1 ~ fit ■ PM ~ *x) ~ 7[r 2 (t») ~ Mtx)] ■ P] 


( 2 . 122 ) 


(2.123) 


This is the expression for the geometric delay that would be observed in the geocentric (S') frame in 
terms of the geometric delay and station positions measured in the Solar System barycentric system 

(£)• 

Since our calculation starts with station locations given in the geocentric frame, it is convenient 
to obtain an expression for [ra(ii) — i“i (t 1 ) ] in terms of quantities expressed in the geocentric frame. 
To obtain such an expression consider two events [r'^tj), 1*2 (^ 1 )] that are geometrically separate, but 
simultaneous, in the geocentric frame, and occurring at time t These two events appear in the Solar 
System barycentric frame as: 



rx(ti) = rift) + (7 - lK(ti) • ^ + lPt\ 

(2.124) 

and as: 

P2(«2) = P(,(ti) + (7 - l)ri(t'i) • ^ + lPt\ 

(2.125) 

where 

t2-tx = 7[pi(ti)-ri(ti)] P 

(2.126) 

With these three equations and the expression 



* 2 (^ 2 ) = P a (tx) + - *x] 

(2.127) 


we may obtain the vector ra(ti): 


r 3 (ti) = r^(ti) + (7 - l)*2(ti) ' 




(2.128) 
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This is the position of station #2 at the time t\ as observed in E. From this we obtain: 


r 2 (ti) -ri(ti) = r^(ti) -r'Jti) + (7 - l)[r^(ti) - r'Jt'J] 


00 

0 2 


(2.129) 


As shown in section 2.1 above, all that is needed to obtain t 2 — ti for the case of plane waves are the 
vectors [r 2 (tx) — ri(ti)] and 0 2 - For curved wave fronts we will need to know the individual station 
locations in the barycentric frame as well. These we obtain from (2.124) and (2.128) with t[ set equal 
to zero. Setting t\ = 0 is justified since the origin of time is arbitrary when we are trying to obtain 
time differences. 

In the actual coding of these transformations, the relationship for the transformation of velocities 
is also needed. Taking differentials in (2.117) and (2.118) we have: 


dr = dr' + (7 — l)dr' • + 7 0dt' 

0 


dt = 7 (dt' + dr' • 0) 

Dividing to obtain dr/dt we obtain for station #2 in the E frame: 


02 




7(l+/?' 2 /») 

For station #2 relative to the geocentric origin, we have from (2.62) and (2.63): 

dU 


02 « (IPN—XYt^e 


where 


w B = 7.2921151467 x 10 -6 rad/sec 


(2.130) 

(2.131) 

(2.132) 

(2.133) 

(2.134) 


is the inertial rotation rate of the Earth as specified in Kaplan (1981), p.12. This is not a critical 
number since it is used only for station velocities, or to extrapolate Earth rotation forward for very 
small fractions of a day (t.e., typically less than 1000 seconds). Actually, this expression is a better 

approximation than it might seem from the form since the errors in the approximation, — — = wg, 

dt 

are very nearly offset by the effect of ignoring the time dependence of PN. 

The assumption of rectilinear motion can be shown to result in negligible errors. Using the plane 
wave front approximation (2.2), we can estimate the error St in the calculated delay if we have made 
an error A0 2 in our value for 0 2 : 


St = k • [r 2 (ti) -ri(t x )] 


1-k (0 2 + A0 2 ) 1-k 0 2 


» tA0 2 


(2.135) 


Further, from (2.132) above, 

A0 2 « A0? 2 

since 


(2.136) 


7 « 1 + 10 -8 (2.137) 

For the vector §t 2 in a frame rotating with angular velocity w, the error A0' 2 that accumulates in the 
time interval r due to neglecting the rotation of that frame is 

A(f 2 « 0 2 w t (2.138) 
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Thus for typical Earth-fixed baselines, where r < 0.02 sec, neglect of the curvilinear motion of station 
#2 due to the rotation of the Earth causes an error of < 4 x 10 -14 sec, or 0.0012 cm, in the calculation 
of r. Similarly, neglect of the orbital character of the Earth’s motion causes an error of the order of 
0.00024 cm maximum. 

The position, Re, and velocity, fi Et of the Earth’s center about the center of mass of the Solar 
System are: 

Re = (2.139) 

fis = (2140) 

2 , rr*i 

where the index i indicates the Sun, Moon, and all nine Solar System planets, rrn is the mass of the 
body indexed by t, while R, and are that body’s center-of-mass position and velocity relative to the 
center of the Earth in the celestial frame. In a strict sense, the summation should be over all objects 
in the Solar System. Except for the Earth-Moon system, each planet mass represents not only that 
planet’s mass, but also that of all its satellites. The R* and are obtained from the JPL planetary 
ephemeris (DE200 as of May, 1982) for the J2000.0 frame. 

In this treatment we have neglected the general relativistic effects which arise from the fact that 
the Earth is spinning in a gravitational potential well. We also have neglected the fact that there are 
differences in the local Earth-caused gravitational potential between the two antenna locations. These 
deficiencies of the model have been estimated to cause errors of less than 0.6 cm in the measurement 
of a 6,000 km baseline (Thomas, 1975). 

Working in a frame at rest with respect to the center of mass of the Solar System causes relativistic 
effects due to the motion of the Solar System in a “fixed frame” to be included in the mean position 
of the sources and in their proper motion. The effect of galactic rotation can be easily estimated. In 
the vicinity of the Sun, the period for galactic rotation is approximately 2.2 xlO 8 years. Our distance 
from the center is approximately 10 kpc = 3.086 XlO 22 cm. Thus, our velocity is 

P = \ ^ « 10 -3 (2.141) 

For a source at zero galactic latitude, the maximum change in apparent position (over one half galactic 
rotation) is 

A0 « ; — - « 6 X 10 -6 arcsec/year (2.142) 

period 

Since a 1-arcsecond angle subtends a distance of 30 meters at one Earth radius, neglecting this effect 
is roughly equivalent to 0.015 cm/year on intercontinental baselines. It also is nearly one to two orders 
of magnitude smaller than the effects of typical source structure, though it is systematic. 


2.11 ANTENNA GEOMETRY 

The above work indicates how the time delay model would be calculated for two points fixed 
with respect to the Earth’s crust. In practice, however, an antenna system does not behave as an 
Earth-fixed point. Not only are there instrumental delays in the system, but portions of the antenna 
move relative to the Earth. To the extent that instrumental delays are independent of the antenna 
orientation, they are indistinguishable to the interferometer from clock offsets and secular changes 
in these offsets. If necessary, these instrumental delays can be separated from clock properties by 
a careful calibration of each antenna system. That is a separate problem, treated as a calibration 
correction ( e.g ., Thomas, 1981), and will not be addressed here. 

However, the motions of the antennas relative to the Earth’s surface must be considered since 
they are part of the geometric model. A fairly general antenna pointing system is shown schematically 
in figure 5. The unit vector, 8, to the apparent source position is shown. Usually, a symmetry axis 
AD will point parallel to s. The point A on the figure also represents the end view of an axis which 
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allows rotation in the plane perpendicular to that axis. This axis is offset by some distance H from a 
second rotation axis BE. All points on this second rotation axis are fixed relative to the Earth. 
Consequently, any point along that axis is a candidate for the fiducial point which terminates this 
end of the baseline. The point we actually use is the point P. A plane containing axis A and perpen- 
dicular to BE intersects BE at the point P. This is somewhat an arbitrary choice, one of conceptual 
convenience. 

Consider the plane Q which is perpendicular to the antenna symmetry axis, AD, and contains 
the antenna rotation axis A. For plane wave fronts this is an isophase plane (it coincides with the 
wave front). For curved wave fronts this deviates from an isophase surface by » H 2 /(2H), where R 
is the distance to the source, and H is taken as a typical antenna offset AP. For H « 10 meters, R 
= Rmoon — 60 Re « 3.6 x 10 8 meters, H 2 /(2R) » 1.4 X 10“ 7 meters, and is totally negligible. R 
has to be 5 km, or 10 _3 JIe, before this deviation approaches 1 cm contribution. Consequently, for all 
anticipated applications of radio interferometry using high-gain radio antennas, the curvature of the 
wave front may be neglected in obtaining the effect on the time delay of the antenna orientation. 



Figure 5. A generalised schematic representation of the geometry of a steerable antenna 

Provided the instrumental delay of the antenna system is independent of the antenna orientation, 
the recorded signal is at a constant phase delay, independent of antenna orientation, at any point on 
the Q plane. Since this delay is indistinguishable from a clock offset, it will be totally absorbed by 
that portion of our model. 

The advantage of choosing the Q plane rather than some other plane parallel to it is that the axis 
A is contained in this plane, and the axis A is fixed relative to the BE axis by the antenna structure. If 
l is the length of a line from P perpendicular to the Q plane, the wave front will reach the Earth-fixed 
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point P at a time At = l/c after the wave front passes through the axis A. If To is the model delay 
for a wave front to pass from P on antenna #1 to a similarly defined point on antenna #2, then the 
model for the observed delay should be amended as: 


t = To — (At 2 — Ati) = ro + ( li — h)/c (2.143) 

where the subscripts refer to antennas #1 and #2. 

For the inclusion of this effect in the model, we follow a treatment given by Wade (1970). Define 
a unit vector I along BE, having a positive sense away from the Earth. Further, define a vector, L, 
from P to A. Without much loss of generality in this antenna system, we assume that s, L, and I are 
coplanar. Then: 

L = ±H i X ' iX !! (2.144) 

|IX|SXl1 . 1 

where the plus or minus sign is chosen to give L the direction from P to A. The plus sign is used if, 
when 8 and L are parallel or antiparallel, the antenna comes closer to the source as H increases. Since 

Ix[sxl) = B-i[Is] (2.145) 

1 = bL = ±H>Ji-[b -if (2.146) 

where the sign choice above is carried through. 

Curvature is always a negligible effect in the determination of a ■ L. Likewise, gravitational effects 
are sufficiently constant over a dimension |L| so as to be able to obtain to a very good approximation 
a single Cartesian frame over these dimensions. Consequently, it is somewhat easier to calculate a 
proper time At = l/c in the antenna frame and to include it in the model by adding it to tq , taking 
into account, in principle at least, the time dilation in going from the antenna frame to the frame in 
which To is obtained. 

Thus, if so is the unit vector to the source from the antenna in a frame at rest with respect to the 
Solar System center of mass, perform a Lorentz transformation to obtain s, the apparent source unit 
vector in the Earth-fixed celestial frame. Actually, the antenna does not “look” at the apparent source 
position s, but rather at the position of the source after the ray path has been refracted by an angle e 
in the Earth’s atmosphere. This effect is already included in the tropospheric delay correction (Section 
4); however since the antenna model uses the antenna elevation angle Eq, the correction must be made 
here as well. For the worst case (elevation angle of 6°) at average DSN station altitudes, the deflection 
can be as large as 2 x 10~ 3 radians. Thus, 61 « He « 0.2 cm for H = 10 meters. A model option 
permits modification of So to take atmospheric refraction into account. The large-elevation-angle 
approximation is the inverse tangent law: 

A E = 3.13 x 10 -4 /tan£ 0 (2.147) 

where E is the elevation angle, and A E the change in apparent elevation induced by refraction. This 
model was implemented only for software comparison purposes, since it gives incorrect results at low 
elevation angles. In the notation of Section 4, a single homogeneous spherical layer approximation 
yields the bending correction 

A E = cos” 1 [cos (Eo + a 0 )/(l + Xo)] - <*o (2.148) 

where 

Xo = ( pz dry + Pz w .,/Mooi)/ & 


ao = cos x [(l -t- <?')/(! + a)] 


a = A/R 


a' = [(l + + 2)/ sin 2 E<^ — lj sin 2 Eq 


(2.149) 

(2.150) 

(2.151) 

(2.152) 


elevation, while the 
1 °. 


This formula agrees with ray-tracing results to within 1% at 6° and ml5% at 1‘ 
corresponding comparisons for Eq. (2.147) give »25% at 6° and a factor of 3 at 

Since we are given I in terrestrial coordinates, we first perform the coordinate transformation 
given by Q above: 

I = Qherrctrial (2.153) 
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With this done, obtain At = l/c , as shown in figure 6 for each of the major antenna types. Note 
that for “nearby” sources we also must include parallax (e.g., geographically separate antennas are 
not pointing in the same direction). If Ro is the position of the source as seen from the center of the 
Earth, and r is the position of a station in the same frame, then the position of the source relative to 
that station is 

R=Ro~r (2.154) 

and in (2.146) we make the substitution 



[Ro — r] ■ I 
|Ro-r| 


(2.155) 


l = ± H \J 1 - (S4) 2 
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S 

o 

1 

c. 



§ 
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l = ± H\j 1 - (Si) 2 


Figure 6. Schematic representations of the four major antenna geometries used in VLBI 


In the modeling software is the facility to provide a time-invariant offset vector in local geodetic 
coordinates (east, north, and local geodetic vertical) from this point (antenna location) to a point else- 
where, such as a benchmark on the ground. This is particularly useful in work involving transportable 
antennas which may be placed in slightly different places relative to an Earth-fixed benchmark each 
time a site is reoccupied. In modeling that offset vector, we make the assumption of a locally plane 
Earth-surface tangent to the geoid at the reference benchmark and the assumption that the local 
geodetic vertical for the antenna is parallel to that for the benchmark. With these assumptions there 
is an identity in the adjustments of antenna location with changes derived for the benchmark location. 
The error introduced by these assumptions in a baseline adjustment is approximately A B X ( d/Rs ), 
where A B is the baseline adjustment from a priori, d is the separation of the antenna from the 
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benchmark, and Re is the radius of the Earth. To keep this error smaller than 0.01 cm for baseline 
adjustments of the order of 1 meter, d < 600 meters is required. 

More troublesome is that an error in obtaining the local vertical by an angle 50, when using 
an antenna whose intersection of axes is a distance, H, above the ground, can cause an error of 
HainSS « HS © in measuring the baseline to the benchmark (Allen, 1982). Unless this error is 
already absorbed into the actual measurement of the offset vector, care must be taken in setting up 
the antenna so as to make 50 minimal. For a baseline error <0.1 cm, and an antenna height of 10 
meters, 50 < 20 arcseconds is required. Often plumb bobs are used to locate the antenna position 
relative to a mark on the ground. This mark is, in turn, surveyed to the benchmark. Even the 
difference in geodetic vertical from the vertical defined by the plumb bob may be as large as 1 arc 
minute, thus potentially causing an error of 0.3 cm for antennas of height 10 meters. Consequently, 
great care must be taken in these measurements, particularly if the site is to be repeatedly occupied 
by antennas of different sizes. 

Another physical effect related to antenna structures is the differential feed rotation for circularly 
polarized receivers. Liewer (1985) has calculated the phase shift 6 for various antenna types. It is 
zero for equatorially mounted antennas. For altazimuth mounts, 

tan 6 = cos <f> sin H / (sin <f> cos 5 — cos <f> sin 5 cos H) (2.156) 

with 4> = station latitude, H = hour angle, and 5 = declination of the source. For X-Y mounts, two 
cases are distinguished: orientation N — S or E — W. The respective rotation angles are 

tan(— 6) = sin <f> ain H / (cos </> cos 6 + sin ^simScos/f) (TV — 5) (2.157) 

tan(— 0) = — cos (sin 5 sin //) (E-W) (2.158) 

The effect cancels for group delay data, but can be significant for phase delay data. The effect on 
phase delay is 

t — (02 — Oi)/f (2.159) 

where / is the observing frequency and the phase rotation at station *'. The feed rotation correction 
is now an optional part of the MASTERFIT model. 
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2.12 PARTIAL DERIVATIVES OF DELAY WITH RESPECT TO 
GEOMETRIC MODEL PARAMETERS 

With respect to any given parameter, the calculation of the time-delay model must be at least 
as accurate as the data is sensitive to that parameter. Consequently, such effects as the curvature of 
the wave fronts were considered. However, such detail is not necessary for determining the derivatives 
with respect to the relevant model parameters. Here, the plane wave approximation is sufficient. 
Iteration on the “solve-for” parameters and the rapid convergence of an expansion of the time delay 
in the relevant parameters about some a priori point permit this simplification. 

In this plane wave approximation we wish to obtain the parameter derivatives with respect to: 

1. the nominal baseline components (actually, station locations), 

2. the parameters of the whole Earth orientation matrix Q described above, 

3. the solid-Earth tidal parameters, 

4. the parameters of source location (right ascension and declination), 

5. the antenna axis offsets, 

6. the constant, ippni in the retardation of the light ray due to gravitational effects. 

The expressions for these derivatives are considerably simplified if tensor notation, with the 
Einstein summation convention, is employed. Before proceeding any further, we make the following 
definitions for this section: 

r = time delay modeled in the geocentric frame, 

r, = this same time delay, but modeled in the Solar System center of mass frame, 

s — source unit vector (in the celestial system at rest with respect to the Solar System center of 

mass), 

0 = velocity of the geocentric frame as measured in the Solar System center of mass frame 
(remember, all distances are measured in time; thus, this quantity is dimensionless), 

02 = velocity of station #2 in Solar System center of mass frame, 

P=l + S-02. This is a factor « 1.0001, which arises from the motion of station #2 during the 
passage of the wave front from station #1 to #2, 

Q = matrix which transforms from the terrestrial system to the celestial system, 

Lo = the baseline vector in the terrestrial system, 

L, = this same baseline vector in the celestial system center of mass frame, 

L = this same baseline vector in the celestial system. 

With these definitions (2.123) may be written 



T — 7(1 — 0 0 2 )u -10 L, 

(2.160) 

For plane waves from (2.2): 

_k-[r 2 -n] b-L, _ s L, 

1 — k • 0 2 1 + 8^2 P 

(2.161) 

Thus, 

T = -i[\-0 2i 0i) 9 -^--i0 k L. k 
P 

(2.162) 

For parameters (represented symbolically by rj) associated with L, k only: 



a,” [ 7<1 +7& ] o,‘ 

(2.163) 

Define the vector: 

Vk = - [7(1 - 02 i^)~ + ^0*] 

(2.164) 
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Then 


dr 

dr) 


r 9L, k 

dtj 


For parameters associated with the source position only: 


dr 

dr) 


= -7(1 -ft, A) — 

P 


ds k 

dr) 


Sk dp 
P dr) 


Since 


P = 1 + *ifti 


Define the vector: 
Then, 

For example: 
Then, 

and 

and 


|^ = -7(l-ft.A)^ 

dr) p 


= -7(l-ft,ft) — 

P 


dak _ 8fcft; dsi 
dr) p dr) 


Ski 


*fcfti 


dsi 

dr) 


My = -7(1 - ft, Pi ) 


L.< 


Sa- 


3 «ft, 


da_ 

da 


da 

d6 


_ »/ ds i 

dr) 1 dr/ 


8 = [ cos 6 cos a, cos 6 sin a, sin 6 ] 


= [ — cos 6 sin a, cos 6 cos a, 0 ] = [ Ai, j4 2 , A 3 ] 


[ — sin 6 cos a, — sin 6 sin a, cos <5 ] = [ Fi,F 2 , F$ ] 


dr 

d^ = MiAi 


dr 

dS 


= MiFi 


Or, if we define the matrices: 

and 

then: 


Ax F, 

G = ( A 2 F 2 

A3 F3 

M — ( Mx, Ms, Ms ) 


= MG 

For station location parameters the algebra is somewhat more complex. Since 


dr dr 
d a' dS 


L, =r 2 (ti) -ri(ti) 

= r 2 (t 2 ) -ri(ti) -fi 2 [t 2 -ti] 

= r 2 (t 2 ) - ri(*i) - 7^ 3 [ r 2( t i) “ r'lK)] • fi 

= te(*'x)-r'x(t'x)j + (7 - l)[ri(t'x) - r'(*i)] fit ~ 7 fi 2 fi • te(*'x) 


(2.165) 

(2.166) 

(2.167) 

(2.168) 

(2.169) 

(2.170) 

(2.171) 

(2.172) 

(2.173) 

(2.174) 

( 2 . 175 ) 

(2.176) 

(2.177) 

(2.178) 

(2.179) 
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we have: 


L, = L + (7 — 1 )L Pfi — 7P2P h 


(2.180) 


or in tensor notation 

Define the tensor: 

Then 

Since 



— Sij + 


(7 ~ 
P 2 



(7 ~ l)ft 

P 2 


-iPii 




L, t — EijLj 

Lj = QjkLo k 
L ti = E{jQjkLo k 


(2.181) 

(2.182) 

(2.183) 

(2.184) 

(2.185) 


Thus, 

r = ^iEijQ jk L 0h (2.186) 

For parameters which are involved with station locations expressed in the terrestrial coordinate system: 


d i r [^E ijQjk ) d -^ = B k d -l± (2.187) 

where the vector element 

B k = ViEijQjk (2.188) 

Such parameters are: r° p . (radius off spin axis), A° (longitude), (height above the equator), r, Pi , 
Aj, ii (their respective time rates), hi (vertical Love number), Z* (horizontal Love number), ipi (phase 
lag of maximum tidal amplitude). The subscript refers to station number, t.e., t = 1, 2. Define the 
matrix: 


W — [— Ri, # 2 ,— Ai, A 2 , —Zi,Z 2 , — Ri , Rz, — Ai, A 2 , —Zi,Z 2 , — Vj, V 2 , — Hi, H 2 , — $ 1 , $ 2 ] (2.189) 


where each column contains the partials of the Lq component vectors x, y, z with respect to the 
parameters. For example, for the constant terms in the station coordinates [see Eqs. (2.37) through 
(2.39)]: 

/ dLo, A 


Ri = 


dr %> 

dL 0v 

dr° 

»Pi 


dU 


0 , 


A< = 


V dr%, J 

( 9L 0l \ 
3A? 
dLo, 
<3A? 
dip, 

V 3A? J 



(2.190) 


~ r %i ®in A° 

r° p .co8A? 


(2.191) 
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For the station coordinate rates, 

Ri = (t - t 0 )Ri A , = (t - t 0 )Ai Zi = (t - t 0 )Zi 

From Eqs. (2.44) through (2.47), and relying on Williams (1970): 


(2.193) 


Vi-\inr\ = s(i)v(i)w(i) o 


(2.194) 


/ d6 ix \ 

9k ( 0 
Hi = ^ = 5(t)K(t)W(t) ?2 (.) 


(2.195) 


(Mi*\ ( dgi(») 

difri drpi 

*■= H - SWVWS'W ^ 

3g3(0 

V 9& / V 


(2.196) 


where t = 1 implies station #1, t = 2 implies #2, and 5(1) = —1, while 5(2) = 1. These partials of 
g with respect to rj> are 

^ = ^ ZkTp ' R,[X * yp “ Y * Ip l (2197) 


dgi. _ 3 n,rl |r 


W 2 4- ,;2 


l / pl - [pp ■ R,][x p X. + s/ p r,] - [n*p - *.s / P ] 2 

\[ x l + yl 


dgz» 3/x, Tp 


9^ J? 

Also, define a vector: 


— pa P ^ [-^*Vp ^p] Y z p + S/p^« , — [2r p ■ R» 


(2.198) 

(2.199) 


dr 

dr 

dr 

dr 

dr 

dr 

dr 

dr 

dr 

dr 

dr 

dr 

“ 9r° 

*Pl 

' dr° 9 

U7 «P2 

9A°’ 

W 

9zf ’ 

dz§’ 

dr, Pl ’ 

dr, P2 ’ 

917’ 

9A? 

9ii ’ 

9z 2 

dr 

9r 9r 

dr 

dr 

dr 

dr 

dr 

1 






dbi ’ 96 2 ’ 9hi ’ 9h 2 ’ 9Zi ’ 9Z 2 ’ 9^1 ’ 5^*2 


D = BW 


(2.200) 

( 2 . 201 ) 


i 
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Certain parameters such as UTl, polar motion, precession, and nutation affect Q only. For these 
parameters 


Define a vector: 


Then 




Ki = * k E ki 


%- = K, ( d -p). 

dr/ \ dr > J 


'0* 


( 2 . 202 ) 

(2.203) 

(2.204) 


for parameters which affect only the orientation of the Earth as a whole. 

A number of parameter partials are available for the orientation of the Earth. These are for 
UTl, X pole, Y pole, and nutation, as well as the angular offset and angular rate terms in the Earth 
orientation perturbation matrix Q. From (2.62): 


Q = n PNUXY 


Define the matrix: 


Y' = 



Then, the partial required for the Y polar motion parameter is: 

dQ 


ae< 


= (IPNUXY' 


(2.205) 


(2.206) 


(2.207) 


An analogous technique is used for the X pole angle, and for UTl, working with the matrix partials 

ax au au 

a©! ’ and a(UTi) aH' 

Partial derivatives of the VLBI observables with respect to the nutation amplitudes and tweaks 
appear formidable at first sight, but are relatively easy to evaluate if the calculation is performed in 
an organised fashion. Symbolizing the parameters by r), we need to evaluate the partials of the matrix 
Q with respect to rf. 

dN TT „du\d6ip v „ 

-U + N — j -^-XY (2.208) 


d, 7 nP i K dSrJ, U + N d6xf> ) dr, 

p [™u + n °!L\*Ji 

dri \dSe dSe J dr) 


Since Sc = e — £, the first partial on the rhs of Eq. (2.209) is equal to 
with respect to Sip and Sc are easily obtained from the expression for N in Eq. (2.80): 


dN 

ds 


(2.209) 


The derivatives of N 


dN 

dSrp 


— sin 6rp 
= I — cos£cos6t/i 


cos e cos 6\p 
■ cos S cose sin Srp 


8inecos 8 ip 
cosSsinesin 6 ip 


and 


dN 

dSe 


— sin l cos Sip — sin S cose sin Sip — sin? sins sin Sip 

— sin e sin S ip cos e sin S ip 

— cos S sin ecos Sip + sin £ cos e cos e cos e cos Sip + sin £ sine 

— sin £ sin e cos Sip — cos £ cos e sin £ cos e cos Sip — cos £ sin e 


From Eq. (2.67), the partials of U with respect to r) are 


au_ 

dr) 


— sin H — cos H 
cos H —smH 
0 0 



( 2 . 210 ) 


( 2 . 211 ) 


( 2 . 212 ) 
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and, from Eq. (2.72), 


(2.213) 


= cos el (cos 2 Sip + cos 2 esin 2 Sip) 
dip ' 

dH 

= —sine tan Sip/(1 + cos 2 etan 2 6^) 

dSe 


(2.214) 


give the required U -dependent terms in Eqs. (2.208) and (2.209). 

Partials of Sip and Se with respect to the parameters Aij and are obtained immediately from 


Eqs. (2.91) and (2.92). For the “free nutations”, 





dSip . 

= sm UfT, 

oAoo 

dSe 

dBoo 

= cos u>fT 

(2.215) 


dSip „ 

— — = T sinw/T, 
oA io 

dSe 

dBio 

= T cos u)fT 

(2.216) 


dSip 

= COSO JfT, 

OA 20 

dSe 

dB^o 

= sin cj/T 

(2.217) 


^±=Tcosw f T, 

0 A 30 

dSe 

dB 30 

= T sinw/T 

(2.218) 

and for the Wahr series terms (j = 1 to 106): 




dSip 

dAoj 

5 

- sin [r 1 

«=x 

dSe 

dBoj 

6 

= cos kfiOiiT)] 

*=1 

(2.219) 

dSip 

dA u 

5 

= Tsin[^A: J< Q!j( r )] > 

1=1 

dSe 

dB xi 

= T cos [y: ky,Q,(r) 
»=i 

(2.220) 

a s ip 

3-A 2 y 

— cos kjiOt* (T)] , 

<=1 

dSe 

dBn 

6 

= sin[^fcy,a,(T)] 
1=1 

(2.221) 

dSip 

3A 3} 

5 

= T cos [y: kjiOti (T)] , 
«=1 

dSe 

dBzj 

6 

= T sin kjiOtj(T) 

*= 1 

(2.222) 


Finally, the partials of the nutation matrix with respect to the “tweaks” A ip and Ae are obtained 
J dN dN 

by making the replacements (2.93) and (2.94) in N. and are then seen to be identical to 

Eqs. (2.210) and (2.211) with the same replacements for Sip and Se. If the a priori tweaks are sero, 
the partials are exactly equal to the expressions (2.210) and (2.211). 

For the parameters in the perturbation matrix, fl, from (2.112): 


an 

360*0 

an 

ase x 


ooo 
= I o o -1 
o 1 

'o o ( 

= I o o -t 
o t i 


(2.223) 


(2.224) 


where t is the number of years from the reference epoch (e.g., J2000). Then, by substituting these 
matrices for n in (2.112), we obtain the appropriate partials of Q for perturbations about the x axis. 
By analogy, the perturbation parameters about the y and s axes may also readily be obtained. 

If we seek the partials of parameters that affect only the “add-on* terms in r = ro + Ar, then 
from (2.123) we have: 


(2.225) 
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for terms which were “added on” in the Solar System barycenter. An example is gravitational bending: 


= 7(1 -ft,) 

dlPPN (1 + Ippn) 

For terms “added on” in the geocentric frame, then: 

dr dA t 
dr) dr) 

An example is the antenna offset vector. In this case: 


dr 

5 (offset station #2) 



and 


dr 

^(offset station #1) 



where the choice of sign for each station is determined by the choice of sign for that 
model portion. 


(2.226) 


(2.227) 


(2.228) 


(2.229) 
station in the 
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SECTION 3 


CLOCK MODEL 


The frequency standards (“clocks") at each of the two antennas are normally independent of each 
other. Attempts are made to synchronize them before an experiment by conventional synchronization 
techniques, but these techniques are accurate only to about 1 /xsec in epoch and 10“ 12 in rate. More 
importantly, clocks often exhibit “jumps” and instabilities at a level that would greatly degrade 
interferometer accuracy. To account for these clock effects, an additional “delay” r c is included in the 
model delay, a delay that models the behavior of a station clock as a piecewise quadratic function of 
time throughout an observing session. Usually, however, we use only the linear portion of this model. 
For each station this clock model is given by: 

To = t c i + r c2 (t - t re /) + r c3 (t - t reJ ) 2 / 2 (3.1) 


The term, t re f, may be set by the user as a specific time (Julian date), or by default taken as the 
midpoint of the interval over which the a priori clock parameters, r cl , r c2 , r c 3, apply. 

In addition to the effects of the lack of synchronization of clocks between stations, there are 
other differential instrumental effects which may contribute to the observed delay. In general, it is 
adequate to model these effects as if they were “clocklike” . However, the instrumental effects on delays 
measured using the multifrequency bandwidth synthesis technique (Thomas, 1981) may be different 
from their effect on delays obtained from phase measurements at a single frequency. 

The bandwidth synthesis process obtains delay from the slope of phase versus frequency 

9 1 

(r = — ) across multiple frequency segments spanning the receiver passband. Thus, any instrumental 
dv 

contribution to the measured interferometer phase which is independent of frequency has no effect on 
the delay determined by the bandwidth synthesis technique. However, if delay is obtained directly 

from the phase measurement, <f > , at a given frequency, u, then this derived phase delay (r p d = 
does have that instrumental contribution. 

Because of this difference, it is necessary to augment the “clock” model for phase delay measure- 
ments: 

Tc pi = To + T c4 (t - t re f) + T cb {t - t ref ) 2 f 2 (3.2) 

where r c is the clock model for bandwidth synthesis observations and is defined in (3.1). Since 
the present system measures both bandwidth synthesis delay and phase delay rate, all of the clock 
parameters described above must be used. However, in a “perfectly” calibrated interferometer, r c 4 
— r c6 = 0. This particular model implementation allows simultaneous use of delay rate data derived 
from phase delay with delay data derived by means of the bandwidth synthesis technique. However, 
our particular software implementation currently is inconsistent with the simultaneous use of delay 
derived from bandwidth synthesis and delay obtained from phase delay measurements. 

To model the interferometer delay on a given baseline, a difference of station clock terms is 
formed: 


^" c 2 To station. 1 


(3.3) 


Specification of a reference clock is unnecessary until the parameter adjustment step, and need not 
concern us in the description of the model. 

The partial derivatives of model delay with respect to the set of five parameters (r c i,r c2 ,r c 3,r c4 ,r c 5) 
for each station are so trivial as to need no further explanation. 
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SECTION 4 


TROPOSPHERE MODEL 


In order to reach each antenna, the radio wave front must pass through the Earth’s atmosphere. 
This atmosphere is made up of two components: the neutral atmosphere and the ionosphere. In 
turn, the neutral atmosphere is composed of two major constituents: the dry and the wet. The dry 
portion, primarily oxygen and nitrogen, is very nearly in hydrostatic equilibrium, and its effects can be 
accurately estimated simply by measuring the barometric pressure. Typically, at sea level in the local 
zenith direction, the additional delay that the incoming signal experiences due to the troposphere is 
approximately 2 meters. Except for winds aloft, unusually strong lee waves behind mountains (e.g., 
Owens Valley, California), or very high pressure gradients, an azimuthally symmetric model based on 
measurements of surface barometric pressure is considered adequate. We have not yet investigated 
where this assumption breaks down, though “back-of-the-envelope” calculations indicate that, except 
in the unusual cases above, the error in such an assumption causes less than 1-cm error in the baseline. 

Unfortunately, the wet component of the atmosphere (both water vapor and condensed water in 
the form of clouds) is not so easily modeled. The experimental evidence (Resch, 1983) is that it is 
“clumpy” , and not azimuthally symmetric about the local vertical at a level which can cause many 
centimeters of error in a baseline measurement. Furthermore, because of incomplete mixing, surface 
measurements are very inadequate in estimating this contribution which even at zenith can reach 20 
to 30 cm. Ideally, this tropospheric induced delay should be determined experimentally at each site. 
This is particularly true for short and intermediate (B < 1000 km) baselines, where the elevation 
angles of the two antennas are highly correlated in the observations. For long baselines, both the 
independence of the elevation angles at the two antennas, and the fact that often the mutual visibility 
requirements of VLBI constrain the antennas to look only in certain azimuthal sectors, allow the use 
of the interferometer data itself to estimate the effect of the water vapor as part of the parameter 
estimation process. For this reason, and because state-of-the-art water-vapor measurements are not 
always available, we also have the capability to model the neutral atmosphere at each station as a 
two-component effect, with each component being an azimuthally symmetric function of local geodetic 
elevation angle. 

At each station the delay experienced by the incoming signal due to the troposphere can be 
modeled using a spherical-shell troposphere consisting of a wet component and a dry component: 


Ttrop station i — ^ wet trop “1“ T dry trop 

The total troposphere model for a given baseline is then: 


(4.1) 


fit — Ttrop station 2 * trop station 1 (4.2J 

If E is the apparent geodetic elevation angle of the observed source at station i, we have (dropping 
the subscript t): 

r trop = PZt r y Rdry (E) + PZ wet ^wet (R') (4.3) 

where pz is the additional delay at local zenith due to the presence of the troposphere, and R is an 
elevation angle mapping function. 

For recent geodetic experiments, the observed delay has been corrected for the total tropospheric 
delay at each station, calculated on the basis of surface pressure measurements for the dry component, 
and water-vapor radiometer measurements for the wet component. This correction is recorded in the 
input data stream in such a way that it can be replaced by a new model. In the absence of such 
external calibrations, it was found that modeling the zenith delay as a linear function of time improves 
troposphere modeling considerably. The dry and wet zenith parameters are written as 

PZi.y, = P°z itW + PZi,* ~ to) (4.4) 

where to is a reference time. 
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Since the model is linear in the parameters p° and p, the partial derivatives with respect to zenith 
delays and rates are trivial. They are: 


dr 




= m** < „ 


(4.5) 


and 


dr 


dpz id 


or w 


(t to)f(i)Jti d w 


where /(*) = 1 for station #2, and —1 for station #1. 


4.1 CHAO MAPPING FUNCTION 


(4.6) 


The simplest mapping function is that obtained by C. C. Chao (1974) through analytic fits to ray 
tracing, a function which he claims is accurate to the level of 1% at 6° elevation angle, and becomes 
much more accurate at higher elevation angles. 


Ri = 


1 


sin E + 


Aj 

tan E -f- Bi 


(4.7) 


where 


A drv = 0.00143 
B dry = 0.0445 
A wet = 0.00035 
B wet = 0.017 


The user must specify values for the zenith delays. 

The partial derivatives of delay with respect to the parameters A drv and B dry are: 


dr 

dA dr y 


-f{i)pz iry R d ry/( tan E + B dry ) 


and 

~QO' d = f{'i)PZdryRdryAdry/[t&nE+ B dry ) 
where R dry is the Chao mapping function, and E is the elevation angle. 


4.2 LANYI MAPPING FUNCTION 


(4.8) 

(4.9) 

(4.10) 

(4.11) 


(4.12) 


(4.13) 


Results of recent analyses of intercontinental data indicate that the Chao mapping function 
[Eq. (4.7)] is inadequate. To rectify this situation, two modifications have been made to the 
MASTERFIT code. First, the dry-troposphere mapping parameters A dry and B dry of the Chao 
mapping function R dry have been promoted to the status of estimable parameters. Second, the code 
now permits the use of a new, more accurate analytic mapping function developed by Lanyi (1984). In 
its simplest form, this mapping function employs average values of atmospheric constants. Provision 
is made for specifying surface meteorological data acquired at the time of the VLBI experiments, 
which may override the average values. Using numerical fits to ray-tracing results, Davis et al. (1985) 
have arrived at another new mapping function, designated the CfA mapping function. Preliminary 
comparisons indicate that the Lanyi and CfA functions are in agreement to within the order of 1 to 2 
cm over a wide range of atmospheric conditions down to 6° elevation angles. Finally, an approximate 
partial derivative is obtained with respect to one parameter in the Lanyi mapping function; this 
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permits adjustment even in the absence of surface data. The Lanyi function was made the default 
MASTERFIT troposphere model in early 1986. 

Motivation for and full details of the development of a new tropospheric mapping function are 
given by Lanyi (1984). Here we attempt to give a minimal summary of the final formulas. The 
tropospheric delay is written as: 

Urop — F(E )/ sin E (4.14) 

where 


F(E) = PZiryFdry{E) + PZy,etFwet{E) 

+ \Pz iTV Fbi{E) + 2 pz iry Pz„. t Fb 2 (E) + pl„' t F b3 (E)\/A + p% iry F M (E)/ A 2 (4.15) 

The quantities pz iry and pz ue , have the usual meaning: zenith dry and wet tropospheric delays. A 
is the atmospheric scale height, A = kTo/mg c , with k = Boltzmann’s constant, To = average surface 
temperature, m = mean molecular mass of dry air, and g c = gravitational acceleration at the center of 
gravity of the air column. With the standard values k — 1.38066 X 10“ 16 erg/K, m = 4.8097 X 10 -23 
g, g c = 978.37 erg/g cm, and the average temperature for DSN stations T 0 = 292 K, the scale height 
A = 8567 m. 

The dry, wet, and bending contributions to the delay, F d ry(E), F we t(E), and Fbi,b2,b3,bi{E), are 


expressed in terms of moments of the refractivity as 

F dry (E) = A 10 (£)G(AM 110 ,u) + 3cruM 210 G 3 (Mi 10 ,u)/4 (4.16) 

F wet {E) = Aoi{E)G(\Mioi/Mooi>u)/Mooi (4.17) 

Fbi(E) = [«tG !3 (A/iio, u)/ sin 2 E — Afo 2 oG 3 (Afi 2 o/Mo 20 , u)]/2 tan 2 E (4.18) 

Fb 2 {E) = — MoiiG 3 (Mhi/Moix, u)/2Mooi tan 2 E (4-19) 

Fb 3 (E) = — Moo 2 G 3 (Mio 2 /Moo 2 > u )/2Af 2 01 tan 2 E (4-20) 

F 64 (E) = Mo 3 oG 3 (Mi 3 o/M 0 3 o, «)/2tan 4 E (4.21) 

A misprinted sign in the last of Eqs. (5) of Appendix B of Lanyi (1984) has been corrected in Eq. 

(4.21). Here G(q,u) is a geometric factor given by 

G(q, u) = (1 + gtt) -1 / 2 (4.22) 


with 

u = 2a j tan 2 E (4-23) 

where a is a measure of the curvature of the Earth’s surface A/R with standard value 0.001345. 

The quantities Ai m (E) and Mu m are related to moments of the atmospheric refractivity, and are 
defined below. Aio(E) involves the dry refractivity, while A 0 i (E) is the corresponding wet quantity. 
The Ai m (E) are given by 


Alrn[E} — M<)(m 


10 n 
+ ££ 


n=l fc=0 


(-l)"+ fc (2n-l)!!M„_ M , m 

u 

n 

Xhdxim 

2"ifc!(n- Jfc)! 

1 "f” XuMiim/Molm 


Xfolm 


(4.24) 


with the scale factor A = 3 for E < 10° and A = 1 for E > 10°. Only the two combinations (l, m) = 
(0,1) and (1,0) are needed for the Ai m (E). The moments of the dry and wet refractivities are defined 
as 

CO 

M nij = J dq q n P dry (q)fl et {q) (4.25) 

o 

where f dr y, wet ( 9 ) are the surface-normalized refractivities. Here, n ranges from 0 to 1, t from 0 to 3, 
and j from 0 to 2; not all combinations are needed. Carrying out the integration in Eq. (4.25) for a 
three-section temperature profile gives an expression for the general moment 

M nii /n\ = (1 - e-^‘)/a" +1 + [l - ^ +n+1 (? 1 .9 2 )] f[ 

i=0 

+ e- a *'T> +n+1 (q u q2)/a n+1 (4.26) 
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Here, 


I 


22(91,92) = 1 - (92 — 9i)A* ( 4 -27) 

The quantities qi and 92 are the scale-height normalised inversion and tropopause altitudes, respec- 
tively. For the standard atmospheric model, 91 = 0.1459 and 92 = 1.424. The constants a and b are 
functions of the dry (a = 5.0) and wet (0 = 3.5) model parameters, as well as of the powers of the 
refractivities (t and j) in the moment definitions. Table IV gives the necessary a’s and 6’s. 


Table IV 

Dependence of the Constants a and b 
on Tropospheric Model Parameters 


t 

j 

a 

b 

1 

0 

1 

a-1 

0 

1 

p 

a0 — 2 

2 

0 

2 

2(a-l) 

1 

1 

0+1 

0(a + 1) - 3 

0 

2 

20 

2{a0 - 2) 

3 

0 

3 

CO 

0 

1 

h- 1 


Note that the normalization is such that M010 = 1; this moment has therefore not been explicitly 
written in Eqs. (4.16) through (4.21). 

At present, provision is made for input of four meteorological parameters to override the standard 
(average) values of the Lanyi model. These are: l) the surface temperature To, which determines the 
atmosphere scale height (default value 292 K); 2) the temperature lapse rate W, which determines 
the dry model parameter a (default values W = 6.8165 K/km, a = 5.0); 3) the inversion altitude 
hi, which determines 91 = hi / A (default value hi = 1.25 km); and 4) the tropopause altitude /12, 
which determines 92 = /12/A (default value hz = 12.2 km). A fifth parameter, the surface pressure 
Po, is not used at present. Approximate sensitivity of the tropospheric delay (at 6° elevation) to the 
meteorological parameters is —0.7 cm/K for surface temperature, 2 cm/(K/km) for lapse rate, and 
—2 cm/km for inversion and 0.5 cm/km for tropopause altitude, respectively. 

Partials of the delay with respect to the dry and wet zenith delays are obtained from Eqs. (4.14) 
and (4.15): 


5 = f{i)[F dry (E) + 2pz iry Fbi {E)J A]/ sin E 

GPZdry 

+ [2pz„ME)/A + 3p| iry F M (^)/A 2 ]/sinE (4.28) 

= /(.•)[£.„(£?) + 2 PZiry F b2 (£)/A + 2pz wet Fbz(E)/ A\/ sin £7 (4.29) 

VPZyuet 

In analysis of data for which meteorological parameters are not available, it is convenient to introduce 
an approximation to the mapping function [Eqs. (4.14) and (4.15)] which involves a one-parameter 
estimate. This parameter p accounts for deviations from standard meteorological conditions. The 
tropospheric delay is expressed as 


. w . _ dr 

t = (Pz irv + Pz „« ) / am E + p — 


(4.30) 
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where the partial derivative is 


5r (pz try + Pz w , t )uMn o 

dp G[Mno, «) [l + G(A/no, u)] sin E 

pz w , t u(Muo — Mioi/Mqqi) 

G(Mno, u)G(Mioi/Mqou u) G(Mno , «) + ^(Afioi/Afooi) «)] sin E 


4.3 CFA MAPPING FUNCTION 


Another approach to improved modeling of tropospheric delay was recently published by Davis 
et al. (1985). Analytic fits to ray-tracing results yield the CfA mapping function 


R< = 


(4.32) 


sin E -+- 


tan E ■ 


sin E -f c 


where E is the elevation angle. The three parameters a, b, c are expressed in terms of meteorological 
data as 



a = 0.0002723 [ 1 + 2.642 x 10 -4 po - 6.400 X 10 -4 e o + 1-337 x 10 -2 T o 
- 8.550 X 10~ 2 o - 2.456 X 10"% ] 
b = 0.0004703 [ 1 + 2.832 x 10 _5 p o + 6-799 x 10 _4 e o + 7.563 x 10 -3 T o 

(4.33) 


- 7.390 x 10~ 2 a - 2.961 x 10 _2 h 3 ] 

(4.34) 


c = - 0.0090 

(4.35) 

Here, po is the surface pressure and eo the surface partial water vapor pressure, both measured in 
millibars. The quantities To, a, and hj have the same meaning and units as in Section 4.2. This 
function is one of three optional mapping functions in the MASTERFIT model. In connection with 
testing parameter estimation for the Lanyi function, the partial derivative of delay with respect to 
surface temperature To in the CfA function was also evaluated. It is 

dr 

PZiryRdry [3-641 x 10 6 (sin E + c)[tan E + 6/(sin E + c)] — 3.557 x 10 6 o] 

(4.36) 

d7b 

(sin E -1- c)[tani?+ 6/(sin E + c)\ 2 
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SECTION 5 


IONOSPHERE MODEL 


The second component of the Earth’s atmosphere, the ionosphere, is a layer of plasma at about 
350 km altitude, created primarily by the ultraviolet portion of the sunlight. In the quasi-longitudinal 
approximation (Spitzer, 1967) the refractive index of this medium is 

n= [ i “(^) 2 ( i± ? cose ) T /2 (5i) 

where the plasma frequency, u p , is 

i/p = (pc 2 r 0 /?r) 1/2 » 8.97 x 10 3 p 1/2 (5.2) 


the electron gyrofrequency, i/ g , is 

_ eB 
9 2irmc 


(5.3) 


and 6 is the angle between the magnetic field B and the direction of propagation of the wave front. 
Here p is the number density of the electrons, and ro is the classical electron radius. 

For the Earth’s ionosphere, with p » 10 12 electrons m -3 , i/ p « 8.9 MHz, while for the interplan- 
etary medium with p w 10 7 — 10 s electrons m -3 , v p » 28 — 89 kHz. In the interstellar medium, 
p » 10 6 electrons m -3 , which gives u p » 3 kHz. At typical microwave frequencies used for geodetic 
VLBI (8.4 GHz), v p /v = 10 -3 for the ionosphere, 10~ 6 for the interplanetary medium, and 3 x 10~ 7 
for the interstellar medium. 

The gyrofrequency, u g , for an electron in the « 0.2 gauss field of the Earth is » 0.6 MHz. Thus, 
for the ionosphere, Vg/v » 2 X 10 -4 at S band (2.3 GHz), and i/„/i/ w7x 10 -6 at X band (8.4 GHz). 
For the interstellar medium B « 10~ 6 gauss, while for the interplanetary region B « 10 -4 gauss. 

Relative to vacuum as a reference, the phase delay of a monochromatic signal transiting this 
medium of refractive index, n, is 




m / 2 4 I 2 4 / 3 

/(*) K(*) +i(*) +-] 


di 


(5.4) 


where 


1 - 1/2 


For 8.4 GHz, we may approximate this effect to parts in 10 6 — 10 7 by: 


Apd « v 2 




1± 


(?) 


1 0 


PS 


-7 


-( 


±{/ g \ 
v / 


cos 6 


(5.5) 


(5.6) 


where 


9 



(5.7) 


and where I e is the t otal numbe r of electrons per unit area along the integrated line of sight. If we 
also neglect the term (i/ 0 cos6)/i/, then the expression for A p d becomes simple and independent of 
the geometry of the traversal of the wave front through the ionosphere: 


A pd = -q/ v 2 


(5.8) 


This delay is negative. Thus, a phase advance actually occurs for a monochromatic signal. Since phase 
delay is obtained at a single frequency, delays derived from phase delay ( e.g ., delay rates) experience 
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a delay increment which is negative (the actual delay with the medium present is less than what 
would be measured without the medium). In constrast, delays measured by means of a group delay 


/ \ . 

technique such as bandwidth synthesis (r = — ) experience an additive delay which can be derived 

av 

from (5.8) by differentiating <f> = i^A p< j with respect to frequency: 


A g d = qjv 2 


(5.9) 


Notice that the sign is now positive, though the group delay is of the same magnitude as the phase 
delay advance. For group delay measurements, the measured delay is larger with the medium present 
than without the medium. 

For a typical ionosphere, r « 1 — 20 x 10 -10 sec at local zenith for v = 8.4 GHz. This effect has 
a maximum at approximately 1400 hours local time and a broad minimum during local night. For 
long baselines, the effects at each station are quite different. Thus, the differential effect may be of 
the same order as the maximum. 

For the interplanetary medium and at an observing frequency of 8.4 GHz, a single ray path 
experiences a delay of approximately 6 x 10“ 7 sec in transiting the Solar System. However, the 
differential between the ray paths to the two stations on the Earth is considerably less, since the 
gradient between the two ray paths should also be inversely proportional to the dimensions of the 
plasma region ( e.g ., one astronomical unit as opposed to a few thousand kilometers). The ray path 
from a source at a distance of 1 megaparsec (3 x 10 7 km) experiences an integrated plasma delay of 
approximately 5000 seconds for a frequency of 8.4 GHz. In this case, however, the typical dimension 
is also that much greater, and so the differential effect on two ray paths separated by one Earth radius 
is still not as great as the differential delays caused by the Earth’s ionosphere. 

These plasma effects can best be removed by the technique of observing the sources at two 
frequencies, is\ and v 2 , where ^ 1,2 v p and where ( 1/2 — Vi |/(i^ 2 + i>i) » 1. Then at the two 
frequencies and u 2 we obtain 

r„i = r + q/vl (5.10) 

and 

Tv2 = T + q/i/% (5.11) 

Multiplying each expression by the square of the frequency involved and subtracting, we obtain 


t — ar„2 + f>T„i 


(5.12) 


where 


and 



(5.13) 


(5.14) 


This linear combination of the observables at two frequencies thus removes the charged particle con- 
tribution to the delay. 

For uncorrelated errors in the frequency windows, the overall error in the derived delay can be 
modeled as 


2 _ 2 _2 1 . 

a r = a a T„ t + b 


v 

Ti/l 


(5.15) 


Modeling of other error types is more difficult and will not be treated in this report. Since the values 
of a and b are independent of q, these same coefficients apply both to group delay and to phase delay. 

If we had not neglected the effect of the electron gyrofrequency in the ionosphere, then instead 
of (5.12) above, we would have obtained 


r = OTi ,2 + 1 + 


q u g cos 0 
U 2 Vi{v2 ~ Vi) 


(5.16) 
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where a and 6 are defined as in (5.13) and (5.14), respectively. 

If we express the third term on the right-hand side in units of the contribution of the ionosphere 
at frequency i/ 2 > we obtain 


r = ar „2 4- br^i 4- 


Aprft/2 t'g cos 9 
v\ (t/2 + Vi) 


(5.17) 


For X band A p< j » 1 — 20 X 10“ 10 sec at the senith. When using S band as the other frequency in the 
pair, this third term is 2 X 10 -4 A P dCos 0 « 2 — 40 XlO -13 sec at zenith. In the worst case of high 
ionospheric electron content, and at low elevation angles, this effect could reach 0.1 cm of total error 
in determining the total delay using the simple formula (5.12) above. Notice that the effect becomes 
much more significant at lower frequencies. 

In the software chain used at JPL, the dual-frequency correction is performed prior to the pro- 
cessing step “MASTERFIT”. MASTERFIT does not have the facility to perform this correction. 
However, the process is described here because it is important to understanding the data input to 
MASTERFIT. For millimeter accuracy, or for lower observing frequencies even at centimeter accuracy 
levels, a correction for the gyrofrequency effect is necessary. 

In the absence of the dual-frequency observation capability described above, one can improve the 
model of the interferometer by modeling the ionosphere, using whatever measurements of the total 
electron content are available. The model we have chosen to implement is very simple. Its formalism 
is very similar to that of the troposphere model, except that the ionosphere is modeled as a spherical 
shell for which the bottom is at the height hi, above the geodetic surface of the Earth, and the top of 
the shell is at the height /i 2 , above that same surface (see figure 7). For each station the ionospheric 
delay is modeled as 

t, = kgI e S(E)/i/ 2 (5.18) 


where 


_ O.lcro 
2 IT 


(5.19) 


/„ is the total electron content at zenith (in electrons per meter squared XlO -17 ), and g — 1(— 1) for 
group (phase) delay. E is the apparent geodetic elevation angle of the source, S{E) is a slant range 
factor discussed below, and u is the observing frequency in gigahertz. 

The slant range factor (see figure 7) is 


J R 2 sin 2 E + 2Rh 2 + hi - J R 2 sin 2 E + 2Rh x + h 2 
S{E) = V (5.20) 

n,2 — ft i 

t 

This expression is strictly correct for a spherical Earth of radius R, and a source at apparent elevation 
angle E. The model employed uses this expression and a geoid surface with a local radius of curvature 
at the receiving station of R equal to the distance from the receiving station to the center of the Earth. 
The model also assumes this same value of R can be used at the ionospheric penetration points, e.g.: 


Ri = R + hi (5.21) 

This is not strictly true, but is a very close approximation, particularly compared to the crude nature 
of the total electron content determinations on which the model also depends. The total ionospheric 
contribution on a given baseline is 


3 r *>t»tion 1 (5.22) 

We assume that the ionospheric total electron content, J e , is the sum of two parts, one obtained by 
some external set of measurements such as Faraday rotation, and the other by some specified additive 
constant: 

le ~ le meat 4" ^ e add (5.23) 

These external measurements, in general, are not along directions in the ionosphere coincident with 
the ray paths to the interferometer. Thus, for each antenna, it is necessary to map a measurement 
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made along one ray path to the ray paths used by the interferometer. Many different techniques to 
do this mapping have been suggested and tried; all of them of dubious accuracy. In the light of these 
problems, and in the anticipation that dual-frequency observations will be employed for the most 
accurate interferometric work, we have implemented only a simple hour-angle mapping of the time 
history of the measurements of J e at a given latitude and longitude to the point of interest. In this 
model we allow the user to adjust the “height” , h, of the ionosphere, but require 


hi = h — 35 km 

h 2 = h + 70 km (5.24) 



Figure 7. The geometry of the spherical ionospheric shell used for ionospheric corrections 

Nominally, this “height” is taken to be 350 km. Setting this height to zero causes the program to 
ignore the ionosphere model, as is required if dual-frequency observations have already been used to 
remove the plasma effects. As in the troposphere model, these corrections can also be incorporated 
into, and passed in the input data stream. Then the user is free to accept the passed correction, 
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and use this model as a small alteration of the previously invoked model, or to remove the passed 
corrections. 

The deficiencies of these ionosphere models for single-frequency observations are compounded 
by the lens effect of the solar plasma. In effect, the Solar System is a spherical plasma lens which 
will cause the apparent positions of the radio sources to be shifted from their actual positions by an 
amount which depends on the solar weather and on the Sun- Earth-source angle. Since both the solar 
weather and the Sun-Earth-source angle change throughout the year, very accurate observations over 
the time scale of a year will be virtually impossible. 

Again, this portion of the model is linear in the parameter I e add- Thus, the partial with respect 
to that parameter is 

dr k /(station #) g(data type) S(E) 

- — 


dl. 


e add 


( 5 . 25 ) 


with /( 2) = 1 and /( 1) = — 1. 
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SECTION 0 


MODELING THE PHASE DELAY RATE (FRINGE FREQUENCY) 


The interferometer is capable of producing several data types: group delay, phase delay, and the 
time rate of change of phase delay. Actually, the time rate of change of group delay is also available. 
However, it is not accurate enough to be of significance for geodetic uses. The models discussed above 
are directly applicable either to group delay or to phase delay. However, the model for the time rate 
of change of phase delay (fringe frequency) must be either constructed separately, or its equivalent 
information content obtained by forming the time difference of two phase delay values constructed 
from the delay-rate measurements as shown below. We chose the latter route since then only models of 
delay are needed. The two phase delay values, r pd (t± A), used to represent the delay-rate measurement 
information content are obtained from the expression 


r pd (t ± A) = r m (t ± A) + r r (t) ± T r A 


( 6 . 1 ) 


where r m (t) is the model used in the delay extraction processing step, r r (t) is the residual of the 
observations from that model, and f r is the residual delay rate of the data relative to that model. 
This modeling for the delay extraction step is covered in Thomas (1981), and is done in analysis 
steps prior to and completely separate from the modeling described in this report. The output of 
those previous steps is such that the details of all processing prior to the modeling described here 
are transparent to this step. Only total interferometer delays and differenced total interferometer 
phase delays (divided by the time interval of the difference) are reported to this step. One of the 
requirements of these previous processing steps is that the model delay used be accurate enough to 
provide a residual phase that is a linear function of time over the observation interval required to 
obtain the delay information. A linear fit to this residual phase yields the value of f r , the residual 
delay rate. Using these two values of r pd , obtained by (6.1) above, the quantity, R, is constructed by 
the following algorithm: 

\T pd {t + A) - r pd (t - A)] 


R = 


2A 


( 62 ) 


This value and the group delay measurement, r gd , are the two data types that normally serve as 
the interferometer data input to be explained by the model described in this report. The software, 
however, also has the option to model phase delay, r pd , directly. In the limit A — ► 0, this expression 
for differenced phase delay approaches the instantaneous time rate of change of phase delay (fringe 
frequency) at time t. In practice, A must be large enough to avoid roundoff errors that arise from 
taking small differences of large numbers, but should also be small enough to allow R to be a reasonably 
close approximation to the instantaneous delay rate. A suitable compromise appears to be A n 2 
seconds. Fortunately, A has a wide range of allowed values, and the capability to model interferometer 
performance accurately is relatively insensitive to this choice. 
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SECTION 7 


PHYSICAL CONSTANTS USED 

In the software that has been implemented we have tried to use the constants recommended by 
the IAU project MERIT (Melbourne et al., 1983). Those that have not been defined in the text above, 
but which have an effect on the results that are obtained using the JPL software, are given below: 


Symbol 

Value 

Quantity 

c 

299792.458 

Velocity of light (km/sec) 

*"o 

2.817938 X 10 - 15 

Classical radius of the electron (meters) 

Re 

6378.140 

Equatorial radius of the Earth (km) 

U>e 

7.2921151467 x 10" 6 

Rotation rate of the Earth (rad/sec) 

f 

298.257 

Flattening factor of the geoid 

h 

0.609 

Vertical Love number 

l 

0.0852 

Horizontal Love number 
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SECTION 8 


POSSIBLE IMPROVEMENTS TO THE CURRENT MODEL 


This section lists areas in which the current model can be improved. 

Variations in the rates of the station clocks in the potential well of the Sun: 

We estimate that the delay error is less than 0.6 cm for measurements of intercontinental 
baselines. As the accuracy of observations improves over the next few years, this portion of 
the model will be improved. 

Earth Orientation Models: 

The models we are using for the orientation of the Earth in space are very nearly the best 
available today. However, there are short-period deficiencies that may be as large as 1 to 2 
milliarcseconds, and longer-term deficiencies of the order of 1 milliarcsecond per year (5 cm 
at one Earth radius). VLBI measurements made during the past few years indicate the need 
for revisions of the annual nutation terms and the precession constant that are of this order 
[Eubanks et al. (1985), Herring et al. (1986)]. The 18.6-year term in the IAU nutation series 
may also be in error, and present data spans are just approaching durations long enough to 
separate it from precession. To provide an improved nutation model, we have implemented a 
MASTERFIT option to use the amplitudes of Herring et al. (discussed at the end of Section 
2.7). This will constitute a temporarily better model of the annual and semiannual nutations 
until the IAU series is officially revised. 

Tidal Effects: 

The octupole component of solid Earth tides is neglected in the present model. Payne (1986) 
estimates that it may contribute as much as 0.5 cm to intercontinental baseline lengths. 

Antenna Deformation: 

Gravity loading and temperature variations may cause variations in the position of the refer- 
ence point of a large antenna that are as large as 1 to 2 cm. Liewer ( 1986) presents evidence 
that these effects cause systematic errors and that their dependence on antenna orientation 
and ambient temperature may be modeled. 

Antenna Alignment: 

Hour angle misalignment of the order of 1 arc minute can cause 1 mm delay effects for DSN 
HA-Dec antennas with 7-m axis offsets. 

Subreflector Focusing: 

For DSN 64-m Cassegrain antennas, allowing the subreflector to slew in order to maintain 
focus changes the path delay by » 4 cm over the 6° — 90° elevation range. Simulations 
(Jacobs, 1987) show that this effect is almost entirely absorbed by the clock epoch and local 
station vertical coordinate parameters. For baselines between two 64-m antennas, this causes 
a potential error of up to 8 cm in length. 

Tropospheric Effects: 

For DSN HA-Dec antennas, the altitude of the declination axis can vary by up to 7 meters. 
This altitude variation, when scaled by a mapping function, causes variations of several 
millimeters in the observed delay. This effect may be detectable in phase delay observations 
currently being analyzed. 
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APPENDIX A 


EQUIVALENCE OF TWO METHODS OF CORRECTING FOR 
GRAVITATIONAL BENDING 

The formulation for differential gravitational retardation should, in the limit of large distances from 
the object producing the gravitational well, produce the same “bending” of the light ray as that 
produced by the more commonly available formula for the bending, i.e.: 


0 = 


2(1 + 7ppat)m 
c 2 R 


(A.1) 


Indeed it does, as haw been shown by Thomas in a private communication (1976). 

Assume the geometry in figure 3 where r, is the direction from the center of the gravitational mass 
to the source, now assumed to be at infinity, and R is the distance of closest approach of the light ray 
to the center of this object. Starting with (2.20), and making use of the two approximations: 


r, + Ti ■ r, n - \rf — R*\ « — - 

and ^ 

?i + ?, = RR/ri 

which hold for a source infinitely far away, we obtain: 


A p = - 


(1 + 1ppn)Hv 


In 


1 + 


Ar ■ R( R/rj) 


(1 + 7 ppn)pp 2Ar • R 


R 


(A.2) 

(A. 3) 

(A.4) 


[*7(2ri)] 

If instead of using the above formulation, we altered the apparent direction of the source by an angle 

2(1 + 'Ippn 


e = 


c 2 R 


(A. 5) 


and calculated the resultant change in delay due to this apparent change in source position, we should 
obtain the same change in delay as A p calculated above. 

From figure 3: 

7 = ~r 7 + QR) (1 - © 2 /2) [A. 6) 

[r, + ts±t| 

A. 

where we have made use of the fact that r, • R = 0. Since 0 < 10 , we have to an accuracy of 10~ 5 

arc seconds that 

= r, + 0R (A. 7) 

Thus, ^ 

r 1 = -Ar r \/c = -Ar -9,/c - Ar • R 0/c (A. 8) 

Since 

r = - Ar -9,/c (A. 9) 

the additional delay, A p , due to the bending is: 


_ 2Ar R(1 + 7 ppn)h p 

p C 2 R 


(A.10) 


which is identical to the expression (A. 4) above for the additional delay due to the differential re- 
tardation. Notice that we have glossed over the motion of station #2 in this treatment. However, a 
complete treatment should give essentially the same result to the accuracy of the terms neglected in 
this treatment. 
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APPENDIX B 


MASTERFIT PARAMETERS 


For the convenience of users of MASTERFIT, Table B.I identifies the names of adjustable parameters 
in the code with the notation of this document. Brief definitions and either references to equations 
(in parentheses) or sections (no parentheses) are also given. 


Table B.I 

Glossary of MASTERFIT Parameters 


Parameter 

MASTERFIT name 

Definition 

Reference 

r »P 

RSPINAX aaaaaaaa 

Cylindrical 

(2.37) 

A 

LONGTUD aaaaaaaa 

station 

(2.38) 

z 

POLPROJ aaaaaaaa 

coordinates 

(2.39) 

r»p 

ORSP/DT aaaaaaaa 

Time rates of 

(2.37) 

A 

DLON/DT aaaaaaaa 

change of 

(2.38) 

z 

DPOL/DT aaaaaaaa 

stn. coords. 

(2.39) 

h, l 

*L0VE # aaaaaaaa 

Love numbers 

(2.44) 

i> 

TIDEPHZ aaaaaaaa 

Tide lag 

(2.50) 

1PPN 

GEN REL GAMMA FACTOR 

PPN gamma 

(2.16) 

a 

RIGHT ASCEN.bsssbsssssbs 

Source RA 

(2.171) 

s 

DECLINATION ssssssssssss 

Source dec. 

(2.171) 

©1,2 

t POLE MOTION 

Pole position 

(2.65), (2.66) 

UT 1 - UTC 

UT1 MINUS UTC 

UT1 - UTC 

2.6 

6®x,y,z 

t AXIS TWEAK OFFSET 

Perturbation 

(2.112) 

6&x,v,z 

f AXIS TWEAK RATE 

coefficients 

(2.112) 

Ao j 

NUTATION AMP LTD PS I cjjj 

Nutation 

(2.87) to 

Mi 

NUTATION AMPLTD PSITcjjj 

amplitudes 

(2.92) 

Mj,sj 

NUTATION AMPLTD PSIA 



Bo, 

NUTATION AMPLTD EPS cjjj 



Bii 

NUTATION AMPLTD EPSTcjjj 



B 2},3j 

NUTATION AMPLTD EPSA 




* Vor H 
t X or Y 
| X, Y or Z 

aaaaaaaa station name 

S88S8S8SS888 source name 
c component: S , C for sine, cosine 

j j j Wahr series term number 


i 
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Table B.I cont. 

Glossary of MASTERFIT Parameters 


Parameter 

MASTERFIT name 

Definition 

Reference 

Tel 

C EPOCH aaaaaaaa 

Coefficients 


Tc2 

C RATE aaaaaaaa 

in clock 


r c 3 

DCRAT/DTaaaaaaaa 

model for 


Tc4 

F OFFSETaaaaaaaa 

delay and 


T c b 

F DRIFT aaaaaaaa 

delay rate 

HI 

PZiry 

DRYZTROPaaaaaaaa 

Dry zenith delay 

(4.3) 

PZ„., 

WETZTROPaaaaaaaa 

Wet zenith delay 

(4.3) 

PZ irv 

DDTRP/DTaaaaaaaa 

Zenith delay 

(4.4) 

PZvci 

DWTRP/DTaaaaaaaa 

time rates 

(4.4) 

Adry 

DRYZMAPAaaaaaaaa 

Chao map 

(4.7) to 

Bdry 

DRYZMAPBaaaaaaaa 

parameters 

(4.11) 

P 

DRYMAPSGaaaaaaaa 

Lanyi map 
parameter 

(4.30) 

To 

SURFTEMPaaaaaaaa 

CfA map 

surface temperature 

(4.36) 

/f add 

Z TECADDaaaaaaaa 

Zenith electron 
content 

(5.23) 


aaaaaaaa station name 
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